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Preface 


This  research  project  is  a  part  of  the  ongoing  effort 
here  at  the  Air  Force  Institute  of  Technology  to  develop 
alternate  methods  for  solving  the  two-dimensional  steady 
state  anisotropic  neutral  particle  transport  equation.  The 
thrust  of  this  research  effort  is  towards  an  accurate  and 
cost-effective  solution  of  the  steady  state  transport  of 
neutrons,  gamma  rays  and  high  energy  x-rays  from  a  lov^ 
altitude  nuclear  burst.  This  problem  which  is  modeled  as 
a  point  source  in  a  two-dimensional  cylindrical  (r,z) 
geometry  with  the  air  ground  interface  included,  is  of 
particular  interest  in  the  areas  of  nuclear  weapons  effects 
and  radiation  physics. 

Presently  the  most  widely  used  computational  methods 
for  solving  the  (air-over-ground)  problem  are  Monte  Carlo 
and  discrete  ordinates.  However,  these  methods  have 
severe  limitations  and  computational  problems.  My  research 
plan  was  to  formulate  and  evaluate  a  solution  technique 
which  did  not  have  these  disadvantages.  A  finite  element 
solution  method  which  is  based  on  a  space-angle  synthesis 
flux  expansion  of  bicubic  splines  and  spherical  harmonics 
was  ciiocon.  The  merits  of  this  solution  technique  were 
examined  ana  a  comouter  algorithm  for  the  nu.Tierical  solu¬ 
tion  of  this  problem  was  developed. 
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David  D.  Hardin,  without  whose  direction,  encouragement  and 
many  hours  of  discussion  and  counselling  this  thesis  would 
not  have  been  possible.  I  am  also  grateful  for  the  support, 
advice  and  encouragement  that  was  provided  by  Dr.  J.  Jones 
of  the  Air  Force  Institute  of  Technology  Mathematics 
Department . 

Finally,  I  wish  to  express  my  appreciation  to  my  wife 
and  daughter  for  their  understanding,  patience  and  constant 
support  throughout  this  project.  To  my  wife,  Cynthia,  I 
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this  thesis. 
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Abstract 


A  finite  element  space-angle  synthesis  ■  solution  of  the 
steady  state  anisotropic  Boltzmann  (transport)  equation  in  a 
two-dimensional  cylindrical  geometry  has  been  developed. 
Starting  from  a  variational  principle  the  Bubnov-Ga 1 erkin 
solution  method  was  applied  to  the  second  order  even  parity 
form  of  the  Boltzmann  equation.  A  trial  function  flux  ex¬ 
pansion  in  bicubic  splines  and  spherical  (surface)  harmonics 
was  used.  A  first  scatter  (collision)  source  and  an  exponen¬ 
tially  varying  atmosphere  was  also  incorporated  into  this 
development . 

Finite  element  space-angle  synthesis  (FESAS)  was 
developed  as  an  alternate  solution  approach  and  an  im.- 
provement  in  co:.Tparison  to  the  methods  of  Monte  Carlo  and 
discrete  ordinates.  FESAS  does  not  have  the  inherent 
characteristics  which  have  produced  the  ray  effect  problem 
in  discrete  ordinates.  Also,  FESAS  may  result  in  lox^;er 
compu ta t iona 1  costs,  than  those  of  Monte  Carlo  and  discrete 
ordinates. 

The  second  order  even  parity  form  of  the  Boltzmann 
equation  v/as  derived  and  shown  to  be  symmetric,  positive 
definite  and  self-adjoint.  The  equivalence  of  a  varia¬ 
tional  minimization  principle  and  the  Bubnov-Ga  lerkin 
method  of  woieh^ed  residuals  was  established.  The  finite 
•‘If'ment  space  antrl'"’  synthesis  svstom  of  .'quations  \-i3s 


v  i  i  L 


expanded  and  a  numerical  computer  solution  approach  was  im¬ 
plemented.  A  computer  program  was  written  to  solve  for  the 
trial  function  expansion  (mixing)  coefficients,  and  also  to 
compute  the  particle  flux. 


THE  APPLICATION  OF  FINITE  ELEMENTS  AND 
SPACE-ANGLE  SYNTHESIS  TO  THE  ANISOTROPIC 
STEADY  STATE  BOLTZMANN  (TRANSPORT)  EQUATION 


I  Introduction 


Background 

Tho  Ai r-Over-Cround  Problem.  The  transport  of  neutrons, 
gamma  rays  and  high-energy  x-rays,  away  from  a  low  altitude 
nuclear  explosion  (air-burst),  is  of  special  interest  in 
assessing  the  vulnerability  and  survivability  of  military 
weapon  systems  and  in  makin<g  radiation  exposure  and  dose 
predictions.  This  neutral  particle  transport  problem  in¬ 
creases  in  complexity  because  of  the  exponentially  varying 
air  density  and  the  air-ground  interface.  A  description  of 
neutral  particle  transport  and,  therefore,  the  air-over- 
grotind  problem  is  given  by  the  Boltzmann  transport  equation. 

Numerical  solutions  to  this  problc-m  alroadv  exist. 

The  main  solution  techniques  are  Monte  Carlo  and  discrete 
ordinates.  However,  discrete  ordinates  and  Monte  Carlo  have 
severe  difficulties  and  disadvantages.  'lo  perform  an  accurat 
'  -c  .ai.ur.iLi  n  r‘:c,u;r-'s  ttu'  luo’  o:  'lours  of 

f  .  ■  i  ■  ^  r  t’ L  *-'*  rt .  i  iiu  t  < '  usoci  loss  ciMiiiuiLor  oxc^"* 

o'.ii  i.oii  t:  i.::io  tixin  M'or.co  Carlo,  -lowover,  it  is  sobjoct  to  a 
cempu  La  L  i  ona  1  difficulty  called  rav  -■ffects  (Ref 


Ray  Effects  and  Discrete  Ordinates.  Ray  effects  are  a 
result  of  the  angular  discretization  of  the  particle  flux  in 
the  discrete  ordinate  method.  It  is  not  a  numerical  problem, 
but  originates  in  the  derivation  of  the  discrete  ordinate 
equations.  In  a  physical  sense  these  equations  only  allow 
source  particles  to  travel  in  specific  directions.  However, 
in  most  practical  problems  these  particles  move  in  all 
directions,  .An  in-depth  analysis  of  the  equations  and  the 
nature  and  reasons  for  ray  effects  can  be  found  elsewhere 
(Ref  1:357). 

Ray  effects  produce  non-physical  distortions  of  the 
angular  flux  in  regions  where  there  are  strong  absorbers, 
localized  sources,  or  high  energy  streaming  particles  (Kef 
2:255-268),  These  distortions  in  the  numerical  formulation 
of  the  discrete  ordinate  method  produce  solutions  which  are 
inaccurate.  The  degree  of  these  inaccuracies  is  dependent 
upon  the  specific  problem  and  the  nature  of  the  absorbing 
media  and  sources.  In  the  air-over-ground  problem  this 
effect  v;ill  bo  significant  because  it  is  essentiallv  a 
strear.iini:  particle  problem  with  localized  first  scatter 
sources  , 

A  considerable  amount  of  work  has  already  boon  done 

in  an  att'^mot.  to  eliminate'  rav  effects  from  tlie  S  equations 

n 

ana  the  uiscrete  ordinate  method.  Ray  effects  can  be  miti- 

■ated  b'.'  the  use  oi  a  fine  anon  In  r  mesh  in  the  finite 

Ji  1  ferencinc  schetiie  of  the  S  oquatio:is.  However,  this 

n 

ai’.proach  increases  l!h‘  computa  t  iona  1  time  and  the  already 


high  computer  costs.  Other  approaches  involve  a  spherical 
harmonic-like  formulation  (Ref  2:255-268),  piecewise  bilinear 
finite  element  approximations  (Ref  3:205-217),  and  space- 
angle  synthesis  with  specially  tailored  trial  functions  (Ref 
4:322-343) . 


The  purpose  of  this  research  project  was  to  d  velop  a 
finite  element  solution  to  the  air-over-ground  problem  by 
using  a  space-angle  synthesis  of  bicubic  splines  in  space 
and  spherical  (surface)  harmonics  in  angle.  Specifically, 
a  solution  of  the  monoenergetic  Boltzmann  equation  in  the 
context  of  the  air-over-ground  problem  was  sought.  Working 
from  a  variational  principle  and  using  a  judicious  choice  of 
trial  functions  the  problem  of  ray  effects  may  be  eliminated 
(Ref  3:214),  Also,  this  judicious  trial  function  choice, 
and  a  Bubnov-Galerkin  solution  method  may  be  more  efficient 
and  less  costly  than  Monte  Carlo  or  discrete  ordinates. 

The  steady  state  solution  of  the  Boltzmann  equation 
with  tirst  scatter  sources,  anisotropic  scattering  and  an 
exponentially  changing  atmosphere  is  desired.  The  problem 
is  formulated  from  a  variational  principle  and  in  a  two- 
dimensional  cylindrical  (r,z)  spatial  geometry  with  the  air- 
ovcr-ground  interface  included.  Fluence  values  as  a 
function  of  two  spatial  (r,z)  and  two  angular  (p,;.;)  variables 
are  sought.  This  is  a  four  dimensional  problem.  (-’ir;ally, 
a  numerical  solution  algorithm  and  the  computer  implemen¬ 
tation  (■'f  this  problem  is  required. 
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Assumptions 

Thcro  are  two  basic  assumptions  which  are  made  in  the 
formulation  of  this  prf'blem.  A  time-independent  (steady 
state)  solution  and  axial  symmetry  is  assvimed.  Because  of 
the  exponentially  chanyiti;-:  air  density  the  air-over-ground 
problem  is  four  dimensional  with  a  spatially  dependent 
(r,z)  solution.  An  assumption  of  axial  symmetry  is  made 
n’ossible  by  inner  inn  the  curvature  of  the  earth.  Ivitiiin 
the  procLem  domain  oi  most  practical  problems  the  curvature 
of  the  earth  is  small  and  can  therefore  be  ignored.  Figure 
1  shows  the  spatial  cylindrical  problem  geometry. 

rile  flux  from  an  air  burse  is  non-zero  :"or  a  fraction 
of  a  second  (microseconds).  Therefore,  particle  fluence 
( number /area )  and  not  flux  is  the  more  iseful  qviantity. 

.A  steady  state  formulation  of  the  air-over-ground  problem 
is  obtained  by  integrating  the  time  dependence  out  of  the 
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Boltzmann  equation.  This  integration  which  is  carried  out 
over  time  limits  when  the  flux  is  zero  produces  a  time 
integrated  or  fluence  equation. 

Development 

In  the  next  chapter  of  this  report  the  problem  equa¬ 
tions  are  presented.  A  finite  element  formulation  of  the 
air-over-ground  problem  is  developed  in  Chapter  III.  The 
Bubnov-Galerkin  method  of  v/eighted  residuals  is  incorpora¬ 
ted  into  this  development.  In  Chapter  IV  a  space-angle 
synthesis  of  bicubic  splines  and  spherical  harmonics  is 
performed.  An  interpolation  of  the  source  terms  is  also 
outlined  in  Chapter  IV.  A  computer  implementation  of  the 
problem  solution  is  examined  in  Chapter  V.  Finally,  con¬ 
clusions  and  recommendations  are  presented  in  Chapter  VI. 


II  The  Problem  Equation 


The  application  of  finite  elements  and  a  variational 
principle  to  the  air-over-ground  problem  and  the  mono- 
energetic  steady  state  Boltzmann  equation  is  not  a  new 
concept  (Ref  5).  As  in  the  work  by  Wheaton  (Ref  5)  only 
the  monoenergetic  problem  will  be  considered.  It  is 
assumed  that  energy  dependence  can  be  easily  incorporated 
into  this  treatment  by  the  use  of  standard  multigroup 
methods.  The  air-over-ground  problem  which  is  in  effect 
the  steady  state  transport  of  neutral  particles  can  be 
described  by  the  Boltzmann  (transport)  equation  and 
appropriate  boundary  conditions  as  follows: 


il*Vd(r,r.)  +  ot(r)d(r 


aS(r,r.Mr)f(r,f2')dn  ' 

+  S(?,r.) 


(1) 


This  is  the  one  speed  monoenergetic  Boltzmann  equation  in 
general  geometry  v\;here 


r  =  the  spatial  position  vector, 

=  a  unit  direction  or  velocity  vector, 
v'  =  gradient  ooerator, 

4>  =  angular  particle  fluence  in  particles/ 
unLt  area/s teradian , 

G(-(r)  =  total  macroscopic  interaction  cross 
section  at  spatial  position  f, 


o r  =  macroscopic  scattering  cross  section. . 

Ihe  proDability  of  a  particle  at  positio.i 
r  and  directio..  '  sea  L  ter  in.  >  into 
■  iirection  a.  It  is  a  tuncLi^a  ot  thr 
acatteriug  angle  '.'■i};'  and  noi  a  function 
of  the  individual  directions  (isotropic 
mi^d  in), 
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(Re£  6:57).  u  is  the  cosine  of  the  angle  formed  by  the 
z-axis  and  the  particle  velocity  vector  Q.  y  is  the  angle 
betii^een  the  planes  formed  by  the  r  vector  and  the  z-axis 
and  that  of  the  r.  vector  and  z-axis. 

The  scattering  properties  of  air  show  a  directional 
dependence  which  is  highly  peaked  in  tlie  forward  direction 
especially  at  the  high  particle  energies  that  exist  in  the 
air-over-ground  problem.  Because  of  this  the  exterior  boundary 
condition  for  this  problem  ’./ill  be  approximated  by  a  vacuum 
boundary  condition: 

=  0  for  on  the  boundary  of  the  problem 

domain  and  f. ‘h  <  0  (d) 

whore  ri  is  tn  ■  outward  unit  normal  to  the  boundary  surface. 

In  physical  terms  this  is  a  non-reentrant  boundary  condition. 

No  particles  are  allowed  to  reenter  the  region  once  they 
leave  it. 

In  two-d imens  iona 1  cylindrical  (r,z)  geometry  there 
is  o’-mmt’trv  in  the  anoie  ■.  This  symmetry  can  bo  written 
as 

=  :(f,%)  for 

k'  r  a  ■: 

(r, 'db; 
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This  symmetry  in  >;  is  shoxv^n  in  Figure  3,  where  the  vectors 


n  and  3  are  perpendicular. 

Note  that  because  of  the  exponentially  varying  air 
density  (in  the  z  direction),  only  azimuthal  symmetry  in  the 


angle  y  is  assured.  There  is  no  symmetry  in  y(Cos  G)  and 
th>'roi'ore  ;(r,')  will  not  be  (’oual  to  The  symmetry 

condition,  Kqs  (3a)  and  (3b),  implies  that 

;'(r,Ujx)  ~  oven  function  in  y  (3c) 

At  tlv'  air  '-round  iriLorfac^.'  :(r,T)  is  oenLinuoiis  except 
at  :  ^  0  (iicf  6;!60),  i.e. 

(r,  )  =  '(r,..)  at  z  =0  and  a  f  0  t4) 

air  ’roiini 

i'e-./ivr  Lin-  'io  r  i '/a  L  i iX  '('r,  )  .-ire  net  conrinuous. 


The  problem  geometry  and  coordinate  systems  as  shown 
in  Figure  2  implies  that  when  p  is  equal  to  zero  the  tngle 
X  must  also  assume  a  value  of  zero.  Therefore,  along  the 
z  axis  (p=0)  the  angle  v  does  not  vary  between  0  and  Itt  , 
This  means  that  there  is  no  x  variation  in  !{i(r,n)  at 
p=0  and  that 

=  0  for  0=0  (5) 

ox 
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Ill  The  Finite  Element  Method 


The  finite  element  method  is  a  mathematical  and 
numerical  technique  for  approximating  the  solution  to 
a  large  class  of  problems.  Initially  i*"  was  developed 
and  used  to  solve  problems  in  stress  analysis  (Ref  7:9). 
Later,  as  the  mathematical  foundation  of  the  method  was 
established  it  gained  widespread  acceptance  and  use  in 
solving  a  larger  class  of  problems. 

Finite  element?  are  an  extension  of  the  Rayleigh- 
Ritz  technique  of  first  recasting  the  problem  in  an  equi 
valent  variational  form  and  then  seeking  a  solution  on 
the  basis  of  an  energy  minimization  principle  (Ref  8:1). 
In  the  Rayleigh-Ritz  method  a  solution  in  the  form  of 
a  linearly  independent  set  of  trial  functions  is  assumed 
These  trial  functions  must  satisfy  the  essential 
boundary  conditions.  The  approximate  or  "best"  solution 
to  the  problem  is  the  linear  combination  of  these  trial 
functions  which  maximizes  (or  minimizes)  the  variational 
principle  ( Cunc t  iona  I  ) .  If  this  linear  combination  of 
functions  is  not  an  extremum  (maximum  or  minimum)  of  the 
functional,  then  the  class  (or  space)  of  trial  functions 
is  expanded  by  the  addition  of  more  functions.  This  ex¬ 
pansion  of  the  trial  function  space  is  continued  until  a 
linear  combination  of  functions  which  is  an  extremum  of 
the  functional  is  obtained. 

The  finite  element  solution  technique  is  similar  to 
that  of  Ra vie igh-R i tz .  The  only  difference  lies  in  the 
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choice  of  trial  functions.  In  finite  elements  the  problem 
domain  is  divided  into  smaller  regions  (grids)  or  elements. 
Each  trial  function  is  usually  associated  with  only  a  few 
elements.  Unlike  Rayleigh-Ritz ,  finite  elements  uses 
trial  functions  which  are  zero  over  parts  of  the  solution 
domain  (a  local  basis).  Also,  the  trial  space  (number  of 
trial  functions)  is  expanded  by  using  more  elements  (mesh 
points)  and  not  by  the  addition  of  a  new  class  of  functions. 
Because  of  these  differences  the  finite  element  method  is 
more  adaptable  towards  a  numerical  (computer)  solution  than 
Rayleigh-Ritz . 

There  are  two  basic  approaches  to  the  finite  element 
formulation  of  a  problem.  One  approach  is  to  find  the 
extremum  of  the  functional  which  originates  from  a  varia¬ 
tional  principle  and  the  calculus  of  variations  (Rayleigh- 
Ritz  with  a  local  basis).  The  other  is  by  the  method  of 
weighted  residuals.  The  method  of  weighted  residuals 
does  not  include  the  use  of  a  variational  principle  or 
the- calculus  of  variations.  In  some  problems  a  variational 
principle  has  not  been  developed  or  may  riot  e>cist,  and 
therefore,  the  Rayleigh-Ritz  approach  cannot  be  used. 
However,  in  such  cases,  the  method  of  weighted  residuals 
can  be  used  to  solve  tlie  problem.  Therefore,  the  method 
■  'u'  ueinhted  r>?sidua]s  can  be  extondo.d  to  a  v/ider  class  of 
'■'alems  tnan  Ra  v  ie  i  t'h-R  L  t  z  . 

rill'  an  tnoci  weighted  residuals  is  another  approach 
r<^r  devi'loping  a  sut  of  (aleebraic)  problem  equations  to 
which  the  i  ini  to  element  method  can  be  applied.  There  ntu? 


three  basic  mathematical  "recipes"  through  which  the  method 
of  weighted  residuals  can  be  developed.  These  are  the  methods 
of  least  squares,  collocation  and  Galerkin.  In  some  problems 
v^?here  a  variational  principle  (functional)  exists  it  can  be 
shown  (Ref  9:735)  that  the  Galerkin  method  of  weighted 
residuals  is  equivaleni  to  Rayleigh-Ri tz .  An  identical  set 
of  matrix  equations  and  therefore  the  same  solution  is 
achieved  by  either  method. 

In  the  following  paragraphs  a  variational  principle  for 
the  air-over-ground  problem  and  the  even  parity  form  of  the 
Boltzmann  equation  is  examined.  A  weak  form  of  the  variational 
principle  and  the  boundary  conditions  are  incorporated  within 
this  development.  Finally,  the  Galerkin  method  of  weighted 
residuals  is  discussed  and  an  equivalence  to  the  variational 
approach  for  the  air-over-ground  problem  is  established. 


Even  and  Odd  Parity  Second  Order  Forms 

In  order  that  a  variational  principle  may  be  used  the 
even  and  odd  parity  forms  of  the  anisotropic  Boltzmann 
ocuntion  and  associated  boundary  conditions  will  bo 
developed.  The  starting  point  oi  this  development  is 
Eqs  (1)  and  (2).  Following  the  derivation  of  Kaplan  and 
Uavis  (Kef  10:166)  and  that  ol  Wheaton  (Ref  3)  the  mono- 
energotic  steady  state  Lransport  i>quatiou  can  bo  written  tn 
terms  of  tne  -i.  vecio  by  changing  to  -ii  in  Eq  (1). 
ihis  gives 


r  ^  :i(r);(r,-:) 


f  • 


3(r.-M  ;  (r  ')d': 


f  u) 


The  even  and  odd  parity  terms  will  now  be  defined  as 
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where 


/s  /\ 

*  > 

A  /S 

o“(r  ) 
a  ( r  ,  n  * ».  ) 


even  parity  fluence 
odd  parity  fluence 
even  parity  source 
odd  parity  source 

even  parity  scattering  cross-section 
odd  parity  scattering  cross-section 


Adding  Eqs  (1)  and  (6),  then  dividing  throughout  by  two 
and  using  the  above  definitions  gives 


+  cc(r)  f  (r  ,n) 


rQr  /N  A  A  A 

^  ^  0  ®  (  r  ,  i'i  *  il  ' )  0  (  r  ,  ir;  ')  d r.  ' 

•  T  A  A 

+  S-(r,...) 


(13) 


Using  the  derivation  of  Wheaton  (Ref  5:8),  which  is  also 
reproduced  in  A[)pendix  A,  the  scattering  kernel  term  in  (13) 
e:>n  vr  ;  t :  -n  ;  . 


(r. 


)  ■  I,  r  ; 
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where  the  even  properties  of  the  even  parity  scattering 
cross-section  and  the  even  parity  fluence  have  been  used 
in  the  derivation  of  Eq  (14).  F.q  (13)  now  becomes 


o.7x(f,3)  +  n^(r)'r(r,fi)  -  I  a- (  f  ,  7  •  5^ ') 'f  (  f  ,  3 ')  d: 


'  ‘«IT 

s"(f .-) 


(15) 


Similarly,  by  substractinq  Fqs  (1)  and  (6)  and  rearranginq 
the  scatterinq  kernel  (See  Appendix  A)  gives 


"i'(r,::)  +  :"^(r)x(f,n)  = 


(f  ,3*3  ')x(f  ,  3  ')d: 


S  (r,  ,) 


(16) 


Eqs  (15)  and  (16)  are  referred  to  by  Kaplan  and  Davis  (Kef  10) 
as  canonical  forms.  The  natural  boundary  condition  Eq  (2) 
c^n  also  be  rewritten  as 


,3) 
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(17) 

and 

(r 

s 
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follows  that 
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and 


-  X(r,^2^ 


(21) 


Also  srid  'i'(r,2)  are  continuous  at  the  air  ground  inter¬ 

face  (when  M  r  0)  but  their  derivatives  are  discontinuous. 

A  further  simplification  is  now  introduced  into  this 
development  by  defining  the  even  and  odd  operators, 
and  as 


GS(r)f(?,P.) 

G''(r)f(r,fi) 


=  aj.(r)f(r,fi)  -  ,9.')d^' 

=  a^(r)f(r,P.)  -  ( r  0  f'r  .P  OdiV 


(22) 

(23) 


where  the  scattering  cross  sections  can  be  expanded  in 
spherical  (surface)  harmonics  (see  Appendix  A)  to  give 


G=>(?)f(?,rn  -  c^(?)f(?,P) 

-  (24) 

where  the  even  parity  Legendre  expansion  cross-section 
is  defined  as 

of(?)  = 

^0  for  li  odd  (25) 

and 


C 

o^(r)  =  I,egendro  expansion  scattering  cross-section 
coeliicients  (Re  I  5:2'2) 


Tilt;  derivatio:; 


(;-)  can  •',>  ;  onnT  in  .Xiinonui;;  . 

I'V.Ul 


It, 


* 


J 


The  odd  parity  scattering  cross  section  can  also  be  expanded 
to  give 


-  (26) 


where  only  odd  expansions  in  I 


are  used  or  is  defined  as 


ao^(r) 


c^(r)  for  £  odd 
0  for  £  even 


(27) 


The  G°  and  G  operators  are  self  adjoint  positive  definite 
(Ref  10:174).  Inserting  these  operators  into  Eqs  (15)  ano 
(16)  they  become 


G°(r)T(r,Q)  =  S°(r,^)  -  n-Vx(r.G) 


and 


=  S’^(r,3.)  -  rcVV(r,3) 


(28) 


(29) 


Solving  for  ^(rjil)  and  x(t,il)  from  Eqs  (28)  and  (29)  pro¬ 
duces 


GS(r) 


)j"^|s-(^,?.)  -  3-vx(r,ro} 

=  ‘j^G^(r;j"*  |s^(r  ,G)  -  i:*VV(r,G)J 


TCr.G')  = 
x(i-‘ 

where  using  the  notation  of  Kaplan  and  Davis 


U  U  1 

K  (r)  =  {G  (r)}  =  inverse  of  the  operator 


k  ’  (  r)  -  ;G  (  r  )  I 


-I 


,-ti ,  . 

G  (r) 


inverse  oL  the  operator 
G-(f) 


(30) 

(31) 


(32) 


(33) 
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A  detailed  mathematical  derivation  of  these  inverse  operators 
is  presented  in  Appendix  A.  They  are  defined  as  (Ref  11:481) 


ry  /s  •-1^1 

K°  (  r )  =  Oj.  ( r )  U 

Ct  ^ 

where  cg(r)  is  defined  by  Eq  (25). 

1 


(34) 


k''(0  =  Ot^(r) 


(35) 


Both  K*^  and  K®  must  be  self  adjoint  positive  definite  since 
they  are  the  inverses  of  and  G°  which  are  both  positive 
definite  and  self  adjoint. 

The  functions  '^jj^(i^)  and  Y^^(r;)  are  the  well  known  nor¬ 
malized  spherical  (surface)  harmonics  (Ref  6:609)  which  are 
defined  as 


(36) 


4ti  ( 2.+m)  ! 


and  Y^^(iq)  being  the  complex  conjugate  of  'Y'j^(f^)  is  defined 
as 


2P,+  1  (  o-ml  I 


4n  ’  (i+m)  ! 

where  a  =  Cosd  and  P‘J'(lO  is  the  associated  Legendre  func- 

;  i  -'!!  s  . 

Lqs  (jo  and  (31)  can  now  be  written  as 


(3  7) 


7(9,13)  =  KS(r)  iS°(r  ,0.)  -  l3*Vx(r,P.)} 


(38) 
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and 


X(r,$^0 


K'^(r)  {S^(r  ,n) 


(39) 


Substituting  Eq  (39)  into  Eq  (28)  produces  the  second  order 
even  parity  form  of  the  Boltzmann  equation 


-P-VK^(r)P-y'l'(r  ,P)  +  G2(r)'i'(?  ,S^)  =  8^(9, P) 

-  P‘VK^(r)S^(r,P)  (40) 


and  inserting  Eq  (39)  into  Eqs  (17)  and  (18)  gives  the  appro¬ 
priate  surface  boundary  conditions  for  Eq  (40)  as 

^(r^.P)  +  K''(r^){S  (r^.P)  -  p.;4'(rg,P)} 

=  0  for  P’n  <  0  (41) 

and 

4'(r  P)  -  K''(r  J{s''(r  P)  -  P-V^r-.P)} 

O  a  ^  b 

=  0  for  fi*n  >  0  (42) 


Eqs  (33)  and  (39)  give  the  second  order  odd  parity  form  of 
the  Boltzmann  equation  which  is 


•>  0/*S/S  /S/N  ll'^  '^/^  «iA/S 

-P-VK^(r)P-Vx(r,P)  +  G'^(r)x(r,P)  =  S^(r,P) 


-K'‘’(r)S-(r,P) 


(43) 


Inserting  Eq  (38)  into  Eqs  (17)  and  (18)  the  surface  boundary 
conditions  for  Eq  (43)  are 


+  KS(r)tsS(?.,P)  -  P-Vx(f  ,u)} 

=  0  for  ;.*n  x  0 

(44) 

(t,.,:) 

.3 

-  K^(?)  iS^(r^,;b  -  P-v^(r^,-:)} 

=  0  for  ^.n  >  0 

(451 
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The  Variational  Principle 


Variational  principles  for  the  monoenergetic  transport 
equation  have  been  found.  These  principles  are  related  by 
a  series  of  transformations  (Ref  10:166).  The  functional 
whose  Euler  equations  are  the  even  and  odd  parity  second  order 
Boltzmann  equations  will  bo  used  in  this  work.  Primarily 
the  even  parity  component  will  be  used  because  it  is  always 
positive,  self-adjoint  and  can  be  integrated  to  give  the 
scalar  flux  or  fluence  (Ref  12:148).  The  odd-parity  flux  which 
can  be  negative  integrates  (over  all  directions)  into  the  net 
current  (Ref  13:12),  Another  "nice"  feature  of  this  even  and 
odd  parity  formulation  is  that  it  produces  a  solution  matrix 
which  is  positive  definite  and  symmetric. 

Defining  the  inner  product  of  two  functions  as 


_n<'.)g(:2)dD 


(46) 


where  '•  means  the  complex  conjugate,  the  functional  which 
corresponds  to  Eqs  (40)  (41)  and  (42)  is  given  as  (Ref  10:169) 


F(u) 


y^{<T.-Vu,K^(.T-Vu)> 
-2<u,s">}dr  -  6[ 


+  <u,cSu>  -  2<n-7u,K^S^> 
1  f:*  h  1  u-d.T}  ds 


(47) 


where  represents  a  surface  integral. 

It  can  be  shown  Ref  (10:169)  that  the  Euler  equation 
(stationary  point)  of  this  functional  is  indeed  the  even  parit:\' 
second  order  lorm  ol  the  Boltzmann  oquatien  and  that  Eqs  (41) 
and  (42)  arc  the  natural  boundary  conditions.  A  similar 
lunctional  for  the  odd  parity  equation  has  also  been  found. 
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In  Appendix  B  the  stationary  point  of  Eq  (47)  is  shown  to 
be  a  minimum  and  the  weak  or  Galerkin  form  is  obtained.  This 
weak  form  is 


,G®'l'>}dr  +  I  n*n  I  n'l'dfids 


J- 


(48) 


where 


n  =  test  or  weight  function 

and 


'i''  =  trial  function 


The  natural  boundary  conditions  Eqs  (41)  and  (42)  are  incorporated 
into  the  weak  form  of  equation  (48).  However,  the  boundary 
condition  at  the  '’round  interface  is  an  essential  condition, 

A  solution  to  the  air-over-ground  problem  can  be  obtained 

from  a  solution  to  Eq  (48)  and  this  essential  boundary  condition. 

The  Galerkin  or  weak  form  of  Eq  (48)  produces  a  solu¬ 
tion  matrix  which  is  positive  definite  and  symmetric.  It 
is  Svmmetric  because  the  test  and  trial  functions  are  the 
same  in  the  Galerkin  solution  method.  The  matrix  is  positive 
definite  because  and  are  positive  definite  and  it  is 
obvious  that  for  the  Galerkin  solution  the  term  <A  • '7n  ,  ( f.  •  v  T ) ' 

is  also  positive  defi’oit’’. 


The  Method  of  '.deiehtod  Residuals 

r'ne  metliou  of  weighted  residuals  is  a  straightforward 
and  simple  prescription  for  solving  a  wide  class  of  problems. 


Unlike  Rayleigh-Ritz  it  does  not  depend  on  a  variational  prin 
ciple  or  the  calculus  of  variations.  However,  in  solving 
certain  types  of  problems  the  use  of  a  variational  principle 
or  the  Galerkin  method  of  weighted  residuals  is  equivalent 
and  they  produce  solutions  which  are  identical.  For  the  air- 
over-ground  problem  and  the  second  order  even  parity  form 
of  the  Boltzmann  equation  it  will  be  shown  that  the  Galerkin 
method  and  the  weak  form,  which  is  given  by  Eq  (43), 
are  equivalent  formulations  of  the  same  problem. 

In  the  method  of  weighted  residuals  an  approximate 
solution  which  is  a  linear  combination  of  trial  functions 
is  assumed  to  exist.  These  trial  functions  are  required  to 
satisfy  the  necessary  boundary  and  continuity  conditions. 

The  approximate  solution,  when  inserted  into  the  problem 
equation,  is  then  required  to  be  an  exact  solution  of 
the  problem  with  respect  to  several  weight  functions  (Ref 
7:106).  The  choice  of  i^eight  functions  determines  whether 
the  method  of  weighted  residuals  is  one  of  collocation, 
least  squares  or  Galerkin.  In  the  Galerkin  method  the 
weight  functions  are  chosen  to  be  the  same  as  the  trial 
functions . 

Applying  the  Galerkin  method  of  weighted  residuals  to 
the  second  order  form  of  the  Boltzmann  equation  is  equivalent 
to  us  111'.':  a  variational  principle.  In  Appendi.x  C  the  weak 
form,  Eq  (43),  is  df’rivv’d  from  a  halerkin  formulation  of 
this  iirobleni.  She  natural  boundary  condition  is  incorporated 
into  this  development  and  an  equation  identical  to  Eq  (48) 


is  produced.  Therefore,  solutions  to  the  second  order  form  of 
the  Boltzmann  equation,  by  using  a  variational  principle  or 
the  Galerkin  method  of  weighted  residuals  are  equivalent. 

This  equivalency  exists  because  the  second  order  even  parity 
operator  of  Eq  (40)  is  positive  definite  and  self-adjoint. 


TV  Space-Angle  Synthesis  of  the 
Even-Parity  Anisotropic  Boltzmann  Equation 

A  space  angle  synthesis  solution  approach  has  already 
been  applied  to  the  air-over-ground  problem.  Roberds  and 
Bridgman  used  "specially  tailored"  angular  trial  functions 
to  solve  the  two  dimensional  steady  state  anisotropic 
Boltzmann  equation  Ref  (4:332).  Space  angle  synthesis  was 
applied  directly  to  Eq  (1);  and  not  to  the  second  order  form 
of  the  Boltzmann  equation.  Miller  e_t  £l.  (Ref  13:12)  have 
applied  phase-space  finite  elements  directly  to  the  isotropic 
second  order  Boltzmann  equation  in  x  -  y  geometry.  Wheaton 
(Ref  5)  has  applied  phase-space  finite  elements  to  the  air- 
over-ground  problem.  However  a  space  angle  synthesis  finite 
element  approach  using  a  flux  expansion  in  bicubic  splines 
and  spherical  harmonics  has  not  been  done. 

Because  of  the  complexity  of  the  air-over-ground  problem 
a  space-angle  synthesis  finite  element  approach  seems  to  be 
justified.  This  solution  technique  might  have  several 
advantages,  some  of  which  are: 

1.  The  elimination  of  ray  effects.; 

2.  The  numerical  advantages  of  finite  elements 
in  combination  with  a  space-angle  synthesis 
approach,  may  be  able  to  better  handle  the 

f our -d imcns ional i ty  of  the  problem; 

3.  Anisotropic  scattering  can  be  easily  handled 
by  a  "wise  and  convenient"  choice  of  angular 
trial  functions  ; 
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4.  The  computational  effort  might  be  reduced, 
without  a  loss  in  accuracy.  It  is  expected 
that  bicubic  splines  will  not  require  a 
fine  spatial  problem  grid  (mesh  size);  and 

5.  The  boundary  conditions  and  a  first  scatter 
source  formulation  can  be  easily  handled. 

The  finite  element  space-angle  synthesis  technique  is 
merely  a  spatial  and  angular  expansion  of  the  even  parity  flux 
by  a  tensor  product  of  polynomial  functions  and  spherical 
harmonics.  In  this  work  a  tensor  product  of  bicubic  poly¬ 
nomial  splines  is  used.  This  expansion  is  the  trial  solu¬ 
tion  which  will  be  used  in  the  finite  element  method.  The 
piecewise  bicubic  spline  expansion  becomes  a  local  basis  in 
the  spatial  (p,2)  variable.  However,  the  spherical  (surface) 
harmonics  which  are  defined  throughout  the  angular  problem 
domain  form  a  global  basis.  Therefore,  this  trial  function 
expansion  has  a  dual  basis  —  a  local  basis  in  space  and  a 
global  basis  in  angle.  This  is  the  finite  element  space 
angle  synthesis  method. 

:'he  starting  point  of  this  development  is  the  weak 
form  of  the  second  order  even  parity  Boltzmann  equation 
and  essential  boundary  conditions.  The  air-over-ground  pro¬ 
blem  is  described  by  the  weak  form  of  liq  (43)  and  the 
svmmeLric  eondition  iiq  (20).  An  essential  bonndarv  connicion, 
at  the  air  ground  interface,  is  that  ':(r,i:)  is  continuous 
at  a  -  0  and  ..  f  0  (see  Fig.  ^).  Howev<'r,  the  derivatives 
of  7(r,d)  are  discontinuous . 


In  the  remainder  of  this  section  the  spatial  and 
angular  trial  functions  will  be  examined.  A  system  of 
coupled  equations  for  the  numerical  solution  of  the  air- 
over-ground  problem  will  be  developed.  A  first  scatter 
(collision)  source  will  be  used  in  this  development. 

The  Trial  Functions 

Because  this  is  a  four  dimensional  problem  with 
anisotropic  scattering  a  numerical  solution  technique  is 
required.  The  finite  element  formulation  of  this  problem 
lends  itself  directly  to  such  a  solution  approach.  However 
an  application  of  four  dimensional  phase-space  finite  elements 
to  Eq  (48)  will  be  very  costly  (computer  costs)  and 
inefficient  (Ref  5:33).  This  is  due  to  the  added  complexity 
of  anisotropic  scattering.  Anisotropic  scattering  increases 
the  bookkeeping  and  computational  difficulties.  A  local 
elemental  basis  in  angle  requires  that  the  scattering  contribu¬ 
tion  to  each  element  must  be  computed  on  an  element  by 
element  basis  for  all  space  and  angle  elements  within  the 
problem  domain.  Therefore  a  four-dimensional  phase-space 
finite  element  formulation  of  this  problem  is  not  a  very 
attractive  or  realistic  approach. 

A  close  examination  of  Eq  (48)  shows  that  the  problem 
operator  is  self-adjoint  positive  definite  symmetric.  This 
allo’.-;.'  -he  use  oi  standard  matrix  iterative  solution  tech¬ 
niques  such  as  Gauss-Seidel  ,  Jacobi  or  Succcuss  iv('  over¬ 
relaxation  (Ref  14:183).  Therefore,  a  numerical  method 
that  includes  a  finite  element  solution  might  be  feasible 


if  the  anisotropic  scattering  contributions  can  be  effectively 
dealt  with.  A  space-angle  synthesis  finite  element  development 
with  a  local  spline  basis  in  space  and  a  global  spherical 
harmonic  basis  in  angle  appears  to  meet  this  requirement. 

A  phase-space  finite  element  problem  formulation  which 
eliminated  the  characteristic  lines  of  the  hyperbolic  discrete 
ordinate  equations  has  been  effective  in  mitigating  ray 
effects  (Ret  3:205).  The  well  known  and  double-Pj^,  equations 
of  nuclear  reactor  physics  have  inherent  elliptic  features 
which  eliminate  ray  effects.  Therefore  space-angle  synthesis 
using  spherical  harmonics  seems  to  represent  an  approach 
which  will  eliminate  ray  effects. 


The  trial  function  expansion  which  will  be  used  in  this 
development  is 


T.Z 


IR  +L  +f. 


i  (  2  ,  P 


L - »  L - 1  Z _ :  Z__.  l,m 


i-i  =  l  ir  =  i  i=0  ir;=-> 


(4S;) 


where  c,z  is  the  spacial  coordinate  dependence  in  cylindrical 
geometry  and  u,X  represents  the  A-  angular  variable  in  Fig.  2. 


'i  (  2  ,  P  ,  y  ,  X ) 


A- 

iz  ,  ir 

j  m 


B,^(2) 


3-,(o) 


even  parity  fluence  at  position  r,z  and  in 
direction  u,x> 

flux  expansion  or  mixing  coefficients, 

cubic  poly norial  spline  in  the  z- 
variablt  (z-spline) , 

cubic  polynomial  c-spline, 

spherical  harmonic  function,  Eq  (51). 
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Figure  4.  Cubic  Spline  With 

Evenly  Spaced  Nodes 


iz,  ir,  i  and  m  are  the  trial  function  summation  indices  and 

IZ  =  total  number  of  z“splines, 

IR  =  total  number  of  p-splines, 

L  =  degree  of  the  spherical  harmonic  expansion. 

The  definition  of  a  cubic  polynomial  splines  (Ref  15:89) 
with  evenly  spaced  nodes  (knots)  is 


(X.^  -x)^-4(x.^^-k)% 


I  -./-X.,  --)  '--(x.— )^-':(X:_  -■:)*,  X:_  4X<X._ 

raxn  ji  tui  -  'rii  iuno^ii'n  is  a:u'wn  in  Fix.  4.  The 


+  isin(mx)} 


(51) 


where 


'"Im 


(^  -  m'  ! 
(!l  +  m) ! 


(52) 


The  trial  function  expansion,  Eq  (49)  can  be  made  to  satisfy 
the  symmetry  of  Eq  (20).  This  is  a  symmetric  condition  in  x 
which  directly  implies  that  the  solution  must  also  be  an  even 
function  in  the  variable  X-  Therefore,  the  angular  trial 
function  expansion  of  Eq  (49)  must  also  be  even  in  x-  Droppine 
the  isin(mx)  term  from  Eq  (51)  and  substituting  into  Eq  (49) 
gives 


i(^,E„P,X) 


IR 

/ 


iz=l  ir=l  t-'O  m=0 


f  j  iz ,  ir  ' 

~~k  i  ,m 


B.  (t)*'s). 


(53) 


where 

ana  by  using  the  or  thonor;:ia  I  properties  of  spherical  harmonics 
the  m  index  begins  at  zero  instead  oi  -I  (see  Appendix  D) . 

The  essential  boundary  condition  at  the  air  ground 
interface  must  also  be  applied  to  Eq  (53).  The  fluence 
c  o  n  L  i  ;;U  :  t  y  rt-qn  ’  I'-jini.-nt  s  can  b-..-  satisiied  by  thib  expansi^ni  in 
bicunic  pelynomial  splines  anu  splicrical  harmonics.  Both 
'■!  th-b  f  line  t )  nr  -  c;.'ii  t  muons  throughoiu  the  proolom 


domn i n  . 


The  cubic  splines  are  also  twice  continuously 


differentiable  but  the  spatial  z-derivative  of  the  solution 
fluence  is  not  continuous  at  the  ground  interface.  However, 
the  z-splines  can  be  modified  to  have  discontinuous  deri¬ 
vatives  at  this  interface  (Ref  16).  A  Double-P.,  or  Yvon's 
method  (Ref  6:163)  can  also  be  used  to  accomodate  the  fluence 
discontinuity  at  p  =  0.  In  this  development  the  air-ground 
interface  will  not  be  included  in  the  problem  domain,  and 
therefore,  this  interface  boundary  condition  will  not  be 
enforced . 

Since  a  Galerkin  method  is  being  used  the  test  or  weight 
functions  are 


'i(z,p,p,X)  = 

Substituting  Eqs  (55)  and  (53)  into  Eq  (48)  produces  Eq 

(58)  where  for  simplicity  the  p,z,y,x  dependencies  are  omitted 

and 


and 


'i  =  B.^(z)B.^(p) 


(56) 


.3  ■  =  d  ■  /' .  ) 


(57) 


1  Kn ' 


1  "ivr/ 


.^d.dsi 


(5X) 
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This  is  a  system  of  coupled  algebraic  equations  where  the 


unknown  quantities  are  the 
can  also  be  written  as 


coefficients . 


Eq  (58) 


NxN  Nxl  Nxl 


^’here 


(59) 


N  =  IZ*IR* (L+1) • (L+2)/2  (60) 

K  =  Coefficient  or  stiffness  matrix,  where 
each  element  is  a  summation  of  terms  in 
the  square  brackets  of  Eq  (58) 

A  =  A-;„  mixing  coefficient  vector 

1  j  A/  m 

S  =  Source  vector  which  is  the  right  hand  side 
of  Eq  (58) 


The  coefficients  of  Eq  (53)  will  be  obtained  from 

a  computer  solution  of  Eq  (59),  These  mixing  coefficients 
can  then  be  substituted  into  Eq  (53)  to  give  the  even  parity 
fluence,  i(z,p,|.,x).  In  a  solution  of  Eq  (59)  the  elements  of 
the  K  matrix  and  S  vector  must  be  computed.  This  computation 
inv'olves  an  evaluation  and  summation  of  the  individual  expanded 
terms  of  Eq  (58).  This  expansion  is  carried  out  in  Appendix 

I.’ 

The  directional  gradient  operator  in  cylindrical  geometry 
is  defined  as  (Ref  6:59) 

/  :  ^  V  l-u''~cosv  2(£±)  “  i  ^{^Vl-u-'  sinx)  +  ud  :■  (61 1 


t'his  is  tile  conservative  term  ot  the  directional  derivative 
in  two  dimensional  (p  ,z)  geometry  with  azimuthal  symmetry. 
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The  and  operators  have  been  defined  in  Eqs  (24)  and 
(35)  of  Chapter  III.  The  scalar  product  of  the  velocity  and 
normal  vectors  is  piven  in  Appendix  E  as 


i>n 


for  the  horizontal  outer  surfaces  (top  or 
bottom)  of  the  problem  cylinder,  and  as 

for  the  vertical  surface  (side)  of  the 
cylinder  (62a) 


where 


'•n_ 


on  the  top  surface,  and 
-p  on  tlie  bottom  surface 


(62b) 


and 


iVn.,  =Vl  cosy  (62c) 

The  normal  unit  vectors  h^  and  h^  are  shown  in  Fig.  ,  Appendix  E. 

Expanding  the  expressions  in  Eq  (5S)  produces  an  integral- 
differential  equation  which  has  twenty-eight  terms  (see  Appendix 
E) .  These  terms,  except  for  the  source  terms,  can  be  easily  se¬ 
parated  into  a  product  of  z,  p,  p  and  y  integrals.  This  is  an 
integral  separation  of  variables  which  is  a  direct  result  of  Eq . 
(53);  where,  it  is  assumed  that  the  solution  can  be  expressed  in 
a  form  where  the  dependent  variables  are  separable.  This  separa¬ 
tion  property  simplifies  the  individual  integrals  which  have  tc 
be  evaluated.  It  allows  for  the  evaluation  of  only  single  inte¬ 
grals  and  not  the  more  complicated  double,  triple  or  quadruple 
integrals.  By  this  separation  of  variables  it  may  be  possible 
to  integrate  most  of  these  single  integrals  ana.  ly  t  ica  1  ly  and 
thus  avoid  a  numerical  integration  process. 

The  Spherical  ilarmonio  Intoe.rals.  liio  usi:  of  a  spherical 
harmonic  angular  trial  function  expansion  was  motivated  b\'  six 
basic  c  on  s  i  do  ivi  i  ivM)  s  ; 


1  •’ 


1.  Because  of  the  fjlobal  nature  of  these  functions  the 
computational  effort  will  be  substantially  reduced, 

2.  Spherical  harmonics  are  well~known  functions  with 
orthonormal  and  symmetric  (odd,  even)  properties. 

3.  The  scattering  cross  sections  are  usually  expanded  in 
spherical  harmonics  (see  Appendix  A). 

4.  Spherical  harmomics  will  produce  a  system  of  equations 
which  are  elliptic  and  invariant  under  continuous 
coordinate  rotations  (Ref  1:362), 

5.  An  analytic  or  closed  form  evaluation  of  the  angular 
integrals  might  be  possible. 

6.  The  even  parity  angular  fluence,  Eq  (53),  can  be 
easily  integrated  to  give  the  total  particle  fluence. 

This  integration  is  carried  out  in  Appendix  I. 

The  term-by-term  expansion  of  Eq  (58)  has  been  partly  carried 
out  in  Appendix  E.  The  resulting  angle  integrals  are  only  de¬ 
pendent  on  the  degree  of  the  spliorical  harmonic  trial  function 
expansion  which  is  used.  They  are  not  dependent  on  the  problem 
parameters  and  therefore  they  can  be  independently  evaluated. 

They  can  be  evaluated  once,  and  thereafter,  used  as  a  part  of 
the  problem  input  data. 

Three  approaches  were  pursued  in  an  attempt  to  evaluate  the 
angle  integrals  which  are  produced  by  this  expansion.  The  first 
approacli  was  to  use  the  orthu'’onal  properties  of  the  associated 
I.egendre  functions  and  the  well-known  properties  of  sines  and 
cosines  to  analytically  evaluate  these  integrals.  However,  be¬ 
cause  of  their  complicated  nature  (see  Appendix  F)  a  closed  1 orm 
aia!’.'!  is  in  t.'- re,  I  i  on  was  nor  oasily  obtaii'.oel  Lor  most  o!  tiit';;i, 

iiu?  second  approach  was  to  use  a  computer  routine  that  can 
•  •'.•a  I  n.  i :  r.iios''  i  n  I  ra  J  s  in  a  svmboiic  or  al'.'.obraic  sense.  Sues 

a  routine'  will  transform  the  integrals  into  algebraic  expressions. 
All'  'O'.'  lac'.'sria  wiiich  w.as  written  I'V  tin'  Massachuso  L  L  s  Ins' 

o  !  io  "  twom,  v/as  n.-'-i  in  this  att-iinj.  ih'wor'^, 
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because  of  the  time  constraint  on  this  research  project  and  the 
need  to  learn  a  new  programming  language  this  approach  was  aban¬ 
doned  . 

Finally,  it  was  decided  to  evaluate  these  integrals  by 
a  numerical  integration  technique.  Because  it  is  possible  to 
separate  the  integration  variables,  the  computational  effort 
can  be  substantially  reduced  by  using  a  single  (one  variable) 
integration  routine.  So  as  to  further  reduce  the  computational 
effort,  Eq  (58)  was  completely  expanded  and  twenty  distinct 
angle  integrals  were  identified.  These  integrals  can  be  found 
in  Appendix  F.  A  ten  point  Newton-Cotes  single  integration 
routine  was  used  to  evaltiate  them.  They  were  evaluated  for 
each  combination  of  the  m,  ,  k  and  n  trial  and  weight  function 
subscripts. 

Bicubic  i’olvno-nial  Splines.  The  spatial  dependence 
of  the  particle  fluence  in  the  air-over-ground  problem  is 
approximated  by  a  product  of  cubic  splines.  The  use  of 
cubic  polynomial  splines  in  the  trial  function  expansion  of 
Eq  (53)  requi.res  the  formation  of  a  tensor  product  space. 

1  h.  i  s  snace  is  t.iaee  u;  of  bicubic  polynomial  splines  which 
are  nrouucLs  oi  ^  and  ic-splines  on  a  rectangular  grid 
(Kol  13:131).  The  exact  shape  of  these  bicubic  splines  are 
ebLait'.ed  iro.i  a  varialienal  principle  or  the  equivalent  method 

I  ail-  '  r-si  h,.i  1  ^  (  .  ,  <• .  ,  a  solution  of  Eq  (58)''. 

p( '  1  \-noi:i .  a  1  splint's  art'  being  used  for  the 

1.  I'h'  '•  lit'  .  .-v  t'l.’i  se  centiaueus  and  form  a 
.  '  ■  1  ■  :  .  -on  '  n.i  t  ;  in'  ’  n  t  "  •  ra  1 

■  ■  :  :  i  -  1  i  -I  . 
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This  reduces  the  number  ot  integrals  which  must  be 
evaluated  and  also  produces  a  sparse  and  banded 
coefficient  matrix. 

2.  A  separation  of  the  p  and  z  integration  variables 
is  possible. 

3.  Third  degree  polynomials  such  as  cubic  splines  have 
a  faster  rate  of  convergence  than  those  of  lower 
degree.  Cubic  splines  are  also  twice  continuously 
differentiable  and  thus  they  are  very  smooth  func-X. 
tions.  For  a  given  problem  mesh  size  splines  will 
produce  a  coefficient  (stiffness)  matrix  that  is 
smaller  but  less  sparse  than  hermites  or  lagrange 
polynomials . 

The  expansion  of  Eq  (58)  with  a  trial  function  of  bicubic 
splines  and  spherical  harmonics  is  carried  out  in  Appendix  E. 

A  further  expansion  of  these  equations  and  a  separation  of  the 
variables  of  integration  produced  seventeen  distinct  c  and  z 
integrals.  These  integrals  which  include  the  space  source 
integrals  are  listed  in  Appendi.x  G.  The  source  integrals  are 
derived  from  an  interpolation  of  the  first  scatter  source  over 
the  entire  spatial  problem  domain. 

The  Source  T.  rms .  A  numerical  solution  to  the  air-over- 
ground  problem  and  the  second  order  Boltzmann  equation  requires 
that  the  source  terms  (right  hand  side  of  Eq  (58))  must  be 
evaluated.  Those  ter.ms  form  the  individual  elements  of  the 
problem  source  vector  in  Eq  (59).  The  even  and  odd  parity 
sources  S'"  amd  S  will  be  defined  as  the  first  scatter  or 
collision  even  and  odd  parity  sources.  The  first  scatter 
source  S(r,.:)  is  the  number  density  of  particles  whici.  leave 
the  burst  point  and  undergo  only  one  coilision  before  being 

,  .  ,  '  .  .  'N  • 

sciittered  into  direction  '  at  position  r.  Streaming  neutrons 
which  leave  the  burst  point  and  do  not  collide  before  reaching 
['o.sition  (r,  :)  arc  not  included  in  the  collision  source. 
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The  use  of  a  first  scatter  source  makes  the  air-over¬ 
ground  problem  more  isotropic.  It  removes  the  strongly 
anisotropic  streaming  particles  from  being  a  part  of  the 
problem  source.  Therefore,  the  solution  fluence  of  Eq  (58) 
will  be  the  scattered  even  parity  fluence  Tg(r,n)  and  not 
the  total  even  parity  fluence  T^Crjil).  The  total  even  parity 
fluence  can  be  defined  as 

=  Tg(r,fi)  +  (63) 

where 

Vj(r,il)  =  streaming  uncollided  particles  at 
position  (r  ,i1)  . 

A  precise  mathematical  definition  of  the  S®  and 
sources  will  now  be  developed.  Also  a  source  interpolation 
procedure  will  be  outlined.  This  source  interpolation  is 
used  in  order  to  simplify  the  source  integrals  of  Appendix 
i:  (E-  40  to  E-46)  . 

The  First  Scatter  Sour  ,  The  even  and  odd  parity 
sources  have  been  defined  as 

s''(r,::)  =  v{S(r,r:)  -  S(r,-A)} 

S-(^,  j  '=  ^'S(?,A)  +  S(r,-A)}  (0) 

IL  S  and  S”  are  first  scatter  source  densities  then  S(r,,:) 
.iiui  S(r,-.)  musi;  also  be  deiiut’d  as  first  scatter  source 
part  ic  li?s/uni  t  volume.  If  a  (losition  (.,•_)  in  the  problem 
domain  is  chosor  then  a  unit  vector  from  the  burst  point 


(0,zb)  can  be  detined  as 


Figure  5.  First  Scatter  (Collision) 
Source  Direction  Vectors 


pCp  +  (z-zb)e^ 
{p^  +  (z-zb)2}^ 


(64) 


Fig,  5  shows  the  direction  vectors  of  this  first  scatter 
(collision)  source,  is  the  direction  that  all  streaming 
(uncollided)  particles  have  at  point  (p,z). 

By  definition  only  particles  which  are  streaming 
radially  outward  from  the  burst  point  can  be  included  in 
the  direct  fluence.  Therefore  the  direct  fluence  at  point 
(p,z)  is  in  the  D'  direction  and 


can  be  written  as 


(65) 


t)j(p  ") 


=  JL,  exp{-f®  a  (z)ds} 

icfr.s"  -'0 


where 


{p2  +  (z-zb)"'}'^ 


(66) 


and  Jds  means  that  the  integration  is  carried  out  along  the 
path  s  (see  Fig.  5).  Also 


Y  =  Total  particle  yield  of  the  nuclear 
explosion  at  the  burst  point. 


at(^)  = 


G|.(0)e“^/®^  for  z  >  0 
a  (ground)  for  z  <  0 


sh  =  atmospheric  scale  height 


(67) 


The  term  J^a^(z)ds  is  the  average  number  of  collisions 
which  a  particle  undergoes  in  traveling  from  the  burst  point 
(o,zb)  to  point  (p,z).  From  Fig.  5  the  distance  s  can  also 
be  written  as 


s  =  (z-zb)/nd  (68) 

and  therefore  by  changing  variables 

ds  =  dz/uil  (69) 


where  ud  is  a  function  of  o  and  z  (but  constant  along  a  path 
lengh  S) 


,z)  -  cos 


(z-zb)/3  =  (z-zD)/  (::-zb)‘}  '  ;70) 


The  integral  term  of  Eq  (65)  can  now  be  written  for  z  >  0  as 

v<°>r: 


(71) 


3S 


and  finally  as 


/  \  “zb/sh  „-z/sh-i 

t(p,z)  =  {e  '  -e  '  }. 


sh 


=  {Ojfzb)  -  a^(z)}*sh 

Ti3’ 


(72) 


From  the  above  derivation  it  follows  that 


t(p,z)  =  { 


{aj.(zb)  -  aj_(z)}*3h/yd  for  z  >  0 
{aj_(zb)  -  Oj.(0)}*sh/yd  -  aj.z/yd  for  z  <  0 


(73) 


also 


'f,(D,2,n")  = 


GXp  (-t(P,z)) 


(74: 


Mote  that  0^(p,z,fi")  is  only  a  function  of  p  and  z. 

The  first  scatter  source  at  (p,z)  and  with  direction 
Q.  are  those  particles  which  undergo  their  first  collision 
at  (p,z)  and  are  scattered  from  direction  ii  '  to  . 
Therefore  the  first  scatter  source  can  now  be  defined  as 


S(p,2,:1)  =  a^(z,3*F.')bj(P,z,^"') 


(75 


where  a  is  not  a  function  of  Q'  but  of  the  scattering  angle 

O  w 

•  'q)  and  z.  From  Fig.  5  and  Fig.  2  a'  is  defined  by 
yd  and  ;<  =  0  i.e.  a'  =  (''',x^)  where  ]i  '  -  yd  and  0. 

By  use  of  the  addition  theorem  it  is  shown  in  Appendix 
H  that  Eq  (75)  can  be  written  as 


s(. 


-z/shV' V 

7  > 

i _ .4-^ 


■iTl 


and  that  the  even  and  odd  parity  first  scatter  sources  are 


L  Z 


S§(P,z.fl)  =  yo=(0)C‘JH-(-l)’''”)Pj„(pd)Pj^(u)ccsmx 

)l=0  nf-=0  (77) 


and 


L 


,  -  s  _  ■  /  sh\  \ 

£=0  nr'>=0 


Oo(0)C^„(l-(-l)'  }Pf^(pd)Pp^(p)cosmx 

(78) 


where  m"  means  that  all  terms  with  a  m  =  0  subscript  must  be 
divided  by  two,  and  cosrax  should  be  interpreted  as  cosine  (mx) 

Source  Interpolation.  Because  of  the  complicated  nature 
of  the  source  expressions  ^  Eqs  (77)  and  (7S),  ai.d  the  need  to 
integrate  the  source  terms  of  Appendix  E,  a  spatial  source 
interpolation  will  be  used.  This  interpolation  which 
simplifies  the  source  integrals  is  necessary  if  a  very 
tedious  (double  or  quadruple)  integration  is  to  be  avoided. 

By  this  interpolation  process  the  source  terms  of  Appendix  E 
can  all  be  separated  into  a  product  of  single  integrals. 

It  is  important  to  note  that  pd  which  is  given  by  Eq  (70) 
is  a  function  of  p  and  z  and  therefore  ,(pd)  is  also  a 
function  of  o  and  z.  Furthermore  <p  ,(P,z,  i^)  of  Eq  (74)  is 
a  fu.nction  of  p  and  z.  Beginning  with  Eqs  (77)  and  (78) 
they  can  be  rewritten  as 


S-^(c 


L  ' 


)F,  (p)cosmx 


(79) 


S^(c,z,n)  = 


o^(0)C:„{l-(-l)  }A.^(p,:r)P.  (y)cosmx 


£=U  m--=Li 


where 


An  fp,z)  =  ;  )P  (;id)e 


-z/sh 


A  spatial  (p,z)  interpolation  of  the  even  and  odd 
parity  sources,  Eqs  (79)  and  (80),  is  therefore  an  inter¬ 
polation  of  A^^(P,2).  In  this  project  these  first  scatter 
second  order  sources  were  interpolated  by  a  combination  of 
piecewise  bi-linear  Lagrange  polynomial  functions.  Speci¬ 
fically,  S°  was  approximated  by  a  tensor  product  of  linear 
Lagrange  polynomials  as  follows 


\  ^  '  O  ' 

S^(pj,2£,i 


>)H.(p)H-(z) 


where 


i=l  j=l 


S^(P',Zj^,L)  =  the  even  parity  source,  Eq  (79) 
^  evaluated  at  the  spacial  nodes 

) 

J 

XZ  =  total  number  of  z-nodes 
NR  =  total  number  of  R-nodes 
P(p)  =  p-linenr  Lagrange  polynM:icl 
L(..:)  ~  ;-l  inear  Lagrange  polync'ml'.l 
jS-;  iin  ;ar  polynomials  arc  delined  as 


n  1  1 1  ,  p  V'  >  —  /  — 


"l-'l 


1 


F  i  g Li re  6  . 


A  Linear  Lagrange  Polynomial  Function. 


The  product  Hj(c)-H£(z)  oi  Eq  (82)  forms  a  tensor  product 
space  on  a  rectangular  grid,  in  the  p,  z  plane  (Ref  15:129) 
Substituting  Eq  (82)  for  in  Eq  (77)  gives 

L 


sS( 


)cm 


NZ  \7l 


X -0  lii'-u 


pm 


Ocosra;/^^ 


i=l  j  =  l  (84) 

■'.imilarly  the  odd  parity  source  can  ulso  be  expressed  as 


S'‘(.  ,z,u)-^  ^  [r^OC^ 

-(-1) 


11 


xm 

Z-TTi^ 


NZ  NR 


P ,  ^(u)cosm' 


an 


'I  ZaJ' i’a)»j(OHi(=)] 


:=i  1=1 


(35) 


ut  ni  ■  ;u thi  -ource  t^rms 
~  ub.;  l  ;  i:  ;i  L  i  on  ninl  a 


■  jo.'arntm  n  ;  •-  :■  r  : r  I  o  s  prouncod  six  distinct 

space  i u t o  r ;i  1  s  .  llu'se  sniico  inlo.irals  are  listed  in  Appendix 


G.  The  source  angle  integrals  have  been  included  in  those 
integrals  which  are  presented  in  Appendix  F. 


V  Computer  Implementation  and  Results 


A  numerical  solution  of  the  air-over-ground  problem 
can  be  obtained  from  a  computer  solution  of  the  system  of 
coupled  algebraic  equations  which  are  represented  by  Eq 

(58) .  These  equations  can  be  written  in  the  matrix  form 
of  Eq  (59)  and  through  a  direct  inversion  or  some  other 
matrix  solution  process  the  j  9  mixing  coefficients 
can  be  found. 

The  computer  solution  to  this  problem  was  accomplished 
by  a  two  step  process.  Each  element  of  the  stiffness  matrix 
and  source  vector  v^as  computed  and  assembled  in  the  matrix 
form  of  Eq  (59).  Then  the  mixing  coefficients  were  computed 
by  the  iterative  matrix  solution  method  of  successive  over¬ 
relaxation.  An  indirect  iterative  matri.x  solution  method 
is  possible  because  Eq  (58)  and  the  stiffness  matrix  of  Eq 

(59)  is  symmetric  positive  definite.  A  computer  program 
which  assembles  the  problem  matrices  and  computes  the 
mixing  coefficient.^  and  particle  fluences  was  written. 

This  pro;rar.i  which  is  written  in  FORTRAN  V  is  listed  in 
.Appendix  J. 

Using  a  ten  point  Nev/ton-Cotes  numerical  integration 
routine,  each  of  the  thirty-seven  integrals  in  Appendix 
F  and  G  were  evaluated.  The  angular  integrals  were  evaluate 
tor  each  .,m,  kn  combination  ot  tlie  trial  anu  weight  tunc  tie. t 
subscrtrjLs.  SitlecLed  prodiict.s  of  theS','  integrals  were  tlum 
used  to  generate  each  ot  the  twenty-eight  terms  (E-i9  to 
E-d6)  of  /\i)[Jc'ndix  E.  FollowinL^  tin:  prescription!  of  .\ppend'.-; 


E  the  first  twenty-one  terms  were  then  added  (or  sub¬ 
tracted)  to  produce  the  elements  of  the  stiffness  matrix. 

The  next  seven  terms  gave  the  elements  of  the  source  vector. 

Writing  the  synthesized  Boltzmann  equation,  Eq  (58),  in 
operator  notation  as 


L( T  •  •  )  n ■  •  =  Sn  . 

''  1  z  ,  1  r  ^  J  2  ,  ]  r  J  2  ,  I  r 

>  ,m  k,n  k,a 

where  T  and  are  defined  by  Eqs  (53)  and  (55)  and 

L(i'-  •  )  ”!  ■  ■  -  Left  hand  side  of  Eq  (58) 

^  12, ir^  J2,jr, 

^  ,  m  k ,  n 

Sn  ■  •  Rivht  hand  side  of  Eq  (58) 

j  a ,  i  r  \  ' 

k ,  n 

the  K  element  of  ttie  stiffness  matrix  is 

p.q 


(B6) 


(37) 


where 


q  =  (ra+1)  +  +  (liiiax  +  l)(Lmax  +  2)‘{(ir-l)  +lR(iz-l)}  (88) 


and 

p  ”  (n-^i)  -  k(k+L)  i  tUiiax  I)(1jkix  +  2)*\(ir-L)  +  JK(iz-l)}  (‘'y) 


The  corresponding  source  vector  element  S,^  is  eiven  bv 
S^,  =  S.B.^(z)B^^(,:.)Q;,^ 
wnore  p  is  given  by  Eq  (89)  and 


(90) 


Lmax  =  degree  of  the  spherical  harmonic  expansion 
JR=IR  =  total  number  of  p-splines 


p  =  the  row  index  of  the  problem  (K)  matrix 
q  =  column  index  of  matrix  K 

B(z),  B(,j)  and  are  defined  in  Chapter  IV  (Eqs  (50)  and 

(54)).  iz,  ir,  jr,  jz,  m,  k,  and  n  are  the  trial  and 
weight  function  expansion  subscripts  where 

f.jk  =  0  to  Lmax 

m  =  Q  to  £ , 
n  =  0  to  k 
iz,jz  =  1  to  17. 
ir,ir  =  1  to  IR 

and 

IZ  =  total  number  of  z-splines 
IR  =  total  number  of  p-splines 

I'sing  the  notation  of  Eqs  (87),  (88)  and  (89)  the  K-problem 
matrix  and  S-source  vector  can  be  easily  assembled  in  the 
following  do  loop. 


DO 

10 

jz  =  1, 

IZ 

DO 

10 

jr  =  1, 

IR 

DO 

10 

k  =  0, 

Lmax 

DO 

10 

n  =  0 , 

k 

S  =  S-.v 
p 

n  • 

1  z  .  1 

DO 

10 

i  z  -  L  , 

lZ 

DO 

iO 

ir  =  1, 

IR 

DO 

10 

4-0, 

Lmax 

DO 

10 

m  =  0 , 
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F. 

p,q 


L(y, 


iz , ir' 


n^. 


.jr 

k,n 


10  continue 


The  assembled  coefficient  (K)  matrix  which  results  from 
Eq  (58)  and  a  bicubic  spline  spherical  harmonic  trial 
function  expansion  is  sparse  and  symmetric.  Because  of  this 
symmetry  and  sparseness  the  coefficient  matrix  can  be  stored 
within  the  computer  by  special  storage  schemes  (Ref  14:70). 
These  sparse  and  symmetric  schemes  will  greatly  reduce  the 
computer  storage  requirements.  As  an  example,  in  a  trial 
function  expansion  where  JR  =  IR  =  8  and  Lmax  =  2  the 
problem  coefficient  matrix  is  a  384  x  384  square  symmetric 
matrix  with  many  elements  which  are  zero.  Therefore  if  this 
entire  matrix  is  stored  within  the  computer  it  will  req.  e 
147456  separate  storage  locations.  This  much  core 
storage  is  already  beyond  the  capacity  of  most  computers.  How¬ 
ever,  by  using  a  sparse  and  symmetric  storage  scheme  this 
matrix  can  bo  reduced  to  one  with  less  than  73728  elements. 

For  large  trial  function  expansions  (JR  =  IR  =  50,  Lmax  =  3), 
aui’cial  auxrliarv  stC)ra'0  ana  solution  tecliniques  will  be 
necessary. 

This  entire  problem  (coefficient  and  source  matrices) 
was  assembled  on  a  CDC  6600  computer  at  the  Air  Force 
Institute  of  Technology.  The  384  x  384  problem  matrix 
is  too  large  to  be  stored  in  core  memory  unless  a  sparse 
symmetric  storage  ’"ode  is  useU.  Because  some  of  the  (p) 
integrals  in  Appendix  G  are  discontinuous  (infinite)  at 
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P  =  0  it  was  necessary  to  use  a  lower  p-integration  limit 
of  1.0  E-8.  Also,  since  the  first  scatter  sources  of  Eqs 
(77)  and  (78)  are  undefined  at  the  burst  point  none  of  the 
problem  nodes  can  be  located  there  (see  Fig,  5). 

Results 

The  computer  routines  which  are  listed  in  Appendix  J 
were  used  to  produce  a  numerical  solution  to  the  air-over¬ 
ground  problem.  These  routines  demonstrate  the  feasibility 
of  using  FESAS  to  produce  a  computer  solution  to  the  two- 
d  Liuens  i  onal  st'>ady  state  anisotropic  Boltzmann  ecuation. 

This  computer  program  has  not  been  fully  developed,  refined  or 
debugged  and  therefore  the  accuracy  of  the  results  has  not 
been  evaluated.  These  results  are  presented  in  an  attempt  to 
further  show  that  FES/vS  is  a  viabli)  numerical  solution 
techn  ique  , 

The  problem  domain  is  a  cylinder  which  sits  on  the  sur¬ 
face  of  the  earth  (see  Fig.  5).  However,  the  air-ground 
interface  is  not  included  in  the  problem  domain  and  there¬ 
fore  all  ground  effocts  are  ignored.  fhe  follov;ing  problem 
parameters  were  usc'.l 

f'eapon  yield  =  l.OE+23  neutrons 

Cylinder  height  =  .4km 

Cvlinder  radius  -  .4km 

Burst  height  -  .12  km 

Total  cross  section  (.'^.(0))  -  15.0  km  ^ 


Table  I 


Legendre  Expansion  Coefficients  which  were  used 
in  a  Numerical  Solution  of  the  Air-Over-Ground  Problem 


Expansion  subscript  1 

Legendre 

expansion  coefficients 

Q 

n 

0 

1 

10.0 

.0 

1 

1  1 

_ _ 

0.0 

2.5 

The  cross-sections  in  Table  I  wore  arbitrarilv  chosen  and  they 
do  not  represent  the  actual  values  for  air  at  sea  level.  A 
relative  convergence  criteria  (.001)  whicli  is  accurate  to  three 
sicniiicant  figures  was  used. 

I  he  pro'iram  execution  times  varied  with  the  degree  of  the 
spherical  harmonic  trial  function  expansion  and  the  problem 
mesh  (grid)  size.  The  entire  problem  matrices  wore  stored 
v/ithin  core  memory  and  by  trial  and  error  it  was  determined 
that  a  relaxation  parameter  of  1.7  gave  the  fastest  convergence 
rate,  Hov;ever,  as  more  trial  functions  were  used  and  the 
system  oi  equations  and  mairices  c.rew  lar'vcr  the  rate  of 
conver ’.enco  decreased.  In  Table  II  the  program  execution  times 
and  the  number  of  iterations  to  convergence  are  listed  for 
.> -ot,  1  ..p;  -nesh  s:::es.  Th<'se  execution  times  and  con- 
■.'■’rn.i’nc’  rates  can  'ne  .mbs  tant  i  a  1 1  v  redne-^d  bv  rev/ritine  or 
■  ie>l  1  t'.-i:’,"  I  li-'  n  r  o'->  r  .''.I'l  routin'''^'  v  vu  pon-d  ;  n  j  . 


FLUENCE^ (NEUTR0NS^KMhh2) 


PARTICLE (NEUTRON)  FLUENCES. 


RflDIUS(KM) 

Figure  7.  Computed  Fluences  as  a  Function  of  Radius  and 
showing  a  variation  with  the  Problem  Spatial 
Mesh  Size. 
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PARTICLE (NEUTRON)  FLUENCES. 


0.2  0.3 

RADIUS [KM) 


Figure  8.  Computed  Pluences  as  a  Function  of  Radius  and 

showing  a  variation  with  the  Degree  of  the  Angular 
Spherical  Harmonic  Trial  Function  Expansion. 
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In  Figures  7  and  8  fluence  values  as  a  function  of  radius  are 
presented.  These  values  are  representative  of  an  altitude  of 
200  metres  and  the  aforementioned  problem  parameters.  Figure 
7  shows  a  variation  of  the  solution  with  spatial  grid  (mesh) 
sizes  whereas  Figure  8  shows  a  variation  with  the  degree  of  the 
spherical  harmonic  trial  function  expansion.  More  detailed 
numerical  results  can  be  found  in  Appendix  K. 


VI  Conclusions  and  Recommendations 


Conclusions 

A  finite  element  space-angle  synthesis  (FESAS)  solution 
of  the  steady  state  anisotropic  Boltzmann  equation  in  two- 
dimensional  cylindrical  geometry  has  been  presented.  In  this 
presentation  a  weak  form  of  the  even  parity  steady  state 
Boltzmann  equation  was  developed.  It  was  shown  that  because 
the  problem  equations  were  positive  definite  and  self-adjoint 
the  Rayleigh-Ri tz  variational  principle  and  the  Bubnov-Galerkin 
method  of  v/eighted  residuals  are  equivalent.  The  problem 
solution  was  formulated  by  using  a  trial  function  expansion  in 
bicubic  polynomial  splines  and  spherical  harmonics.  This  trial 
function  expansion  has  a  dual  basis  —  a  local  basis  in  space 
and  a  global  basis  in  angle. 

This  development  was  specialized  to  the  air-over-ground 
neutral  particle  transport  problem.  It  was  shown  that  a 
finite  element  space-angle  synthesis  solution  is  possible  and 
that  a  first  scatter  interpolation  source  can  be  used.  It 
v.'as  also  sho'wn  that  a  numerical  solutic'n  can  be  achieved  and 
that  this  solution  technique  may  eliminate  ray  efiects  and 
reduce  computational  costs. 

Recomme nda  t i ons 

The  nro  1  i  minary  ro5-,ult3  of  this  st\.idy  have  shown  that 
the  FF.SAS  method  can  produce  a  numerical  solutioTi  to  the 
ateanv  state  Boltc.mann  equation  a:id  t  he  a  i r-ove r-;: round 
problem.  However,  because  of  the  time  constraints  on  this 


research  project  a  complete  development  and  evaluation  of  the 
FESAS  method  was  not  accomplished.  Therefore,  there  are  a 
number  of  recommendations  for  the  further  analysis  and 
evaluation  of  the  FESAS  method,  which  must  be  made.  Some 
of  these  recommendations  are: 

1.  To  enforce  the  boundary  conditions  at  the  air- 
over-ground  interface.  This  can  be  accomplished 
by  a  coalescing  or  stacking  of  the  nodes  (knots) 
of  the  bicubic  splines  and  by  using  a  Double- 

P,,,  approximation  at  this  interface. 

2.  Develop  or  refine  the  computer  algorithm  so  that 
a  comparative  study  can  be  made.  This  study 
should  include  a  comparison  of  the  computational 
costs  and  accuracy  of  FESAS  to  those  of  Monte 
Carlo  and  discrete  ordinates.  Also,  a  determina¬ 
tion  should  be  made  as  to  whether  ray  effects 
have  been  eliminated. 

3.  Obtain,  if  possible,  a  closed  form  solution  to 
the  angle  integrals  in  .-ippendix  F. 

4.  Explore  other  ways  of  handling  the  d i s con t i nu  i  ty 
(at  0  =  0)  of  the  space  integrals. 

5.  1,'so  other  spatial  trial  functions.  Lower  degree 
bi-quadratic  splines,  hermites  and  Lagrange 
polynomials  arc*  possible  candidates.  A  compari¬ 
son  of  the  results  which  arc  obtained  from  the 
use  of  various  trial  functions  can  then  be  made. 

6.  Examine  the  effects  of  an  improved  source  inter¬ 
polation  on  the  sc'lution  accuracy  and  rate  of 
convergence.  .An  improved  source  interpolation 
can  be  achieved  by  the  use  of  a  smalLir  source 
(space)  grid  or  a  higher  degree  Lagrange  or 
Hermito  polynomial  interpolation  function 

7.  Extend  the  use  of  finite  element  space-an'le 
synthesis  to  the  solution  of  energy  dependent 
mnl  t  i-':ronp  problems  . 
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Appendix  A 


Derivation  of  the  Scattering  Kernel,  Inverse  Collision 
Operators  and  the  Even  and  Odd  Parity  Collision  Cross  Sections 


Collision  Cross-Sections 

The  even  and  odd  parity  legendre  expansion  scattering  cross 
sections  a|  and  o'^  originate  in  the  derivation  of  the  second 
order  forms  of  the  Boltzman  equation  (Chapter  III),  A  defini¬ 
tion  of  these  quantities  has  been  given  elsewhere  (Ref  5:29). 
This  definition  is  repeated  below. 

Beginning  with  the  even  and  odd  parity  scattering  cross 
sections 


(A-l) 


(A-2) 


the  usual  computational  practice  (Ref  6:131)  is  to  expand  these 
scattering  cross  sections  in  Le<;,endre  polynomials  as 

Z 

,  'RCA-A.O  (a-3) 

or 


where 

=  macroscopic  scattering  cross  section 

ij^(r)  =  '.ogondre  macroscopic  cross-st:;ctioii 
expansion  coefficient  which  is  a 
fur.ction  of  position  (material) 

L  =  the  degree  legendre  polynomial  expan¬ 
sion  which  is  used 


^A)-2.A±L.~  ECAl-A') 


(A-/f) 


=  Legendre  polynomial  of  degree  Z. 

^•S2'  =  scattering  angle  (u^  =  CosG) 

Inserting  Eqs  (A-3)  and  (A-4)  into  (A-1)  and  (A-2) 


From  the  even  and  odd  properties  of  l.egendre  polynomials 


17:223) 


even  function  if  Z  is  even 
odd  function  if  Z  is  odd 


therefore 


if  Z  is  odd 


s  even 


and 


'2'RC/i-A')  if  i  is 


^~o 


0 


if  Z  is  even 


Eqs  (A-5)  and  (A-6)  can  now  be  rewritten  as 

Z, 

6l(f^2-A')  ^4'-'=)  •  14 (A.  A') 

L—t  4-7r  ^ 


and 


t?y£h 


.2'  -z  / 


(A-5) 

(A-6) 

(Ref 

(A-7) 

(A-8) 

(A-9) 

(A-10) 


(A-11) 


The  even  and  odd  parity  cross-sections  can  also  be 


expanded  in  l.egendre  polynomials  as 


ina 


e=o 


(A-12) 


.  . . 


(A-13) 

Comparing  Eqs  (A-”L3),  (A-12),  (A-11)  and  (A-10)  it  is  apparent 
that 


and 


of(r)  =  \ 


= 


0  j  ('  r )  for  £.  even 
0  for  i  odd 


for  £  even 


<j^(r)  for  £  odd 


(A-14) 


(A-15) 


Therefore  o?(r)  and  o^(r,.Q*i3')  are  even  functions.  Also  cV(r) 


and  a^(r  jSl'il  ■')  are  odd  functions, 
s 


Scattering  Kernel 

In  the  development  of  the  even  and  odd  parity  forms  of  the 
anisotropic  steady  state  Boltzmann  equation  (Chapter  III),  it 
was  necessary  to  express  the  scattering  kernels  in  terms  of 
the  even  and  odd  parity  flux  (fluence)  components.  Following 
the  derivation  of  Wheaton  (Ref  5:11)  the  scattering  terms  can 
be  written  as 


where  the  integrations  are  carried  out  over  all  directions. 
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1 
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Because  is  an  even  function,  meaning  f(x)  =  f(-x), 
it  follows  that 


(A-17) 


and  therefore 


With  the  even  parity  flux  defined  as 


Eq  (A-18)  becomes 


(A-18) 


(A-19) 


and  by  a  similar  derivation 


(A-2G) 


(A-21) 


Inverse  Collision  Op>'rators 

In  deriving  the  second  order  forms  of  the  Boltzmann 
equation  (Chapter  III)  the  G®  and  g'^  operators  were  defined  as 

^  -j (A-22  ) 

wrr 

and 


4?: 


Z.yZ 


'■) 


■^rr 


(A-23^ 


(SI 


where  the  r  dependence  has  been  omitted  in  an  attempt  to 

simplify  the  notation. 

Usin'"  the  addition  theorem  (Ref  6:609) 


1:00  W 


in  (A-12)  and  (A-13)  gives 


Eqs  (A-22)  and  (A-23)  can  now  be  written  as 


(A-24) 


(A-25) 


(A-26) 


(A-27) 
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QYtA) = 

The  inverse  operators  are  defined  as 


/  5 

/c -p 


where  it  is  mean 


t  that  if 


then 


(A-28) 


(A-29; 


(A-3iT, 


multiplying  Eq  (A-30)  by  expanding  the 


operator  gives 


[£i) 


v/r 


/ 

a\  I  //  ^«-\ 


(A-32) 


using  the  orthonormal  properties  of  spherical  harmonics  (Ref 


6:609) 


Ja-fY 


where  Sj,  is  the  F’ronecker  delta  which  is  defined  by 


= 


0  i  f  V. 


1  £  =  k 


Therefore  Eq  (A-32)  can  now  be  written  as 

={^cyt) 


(A-32) 


(A-34) 


(A-35) 


; 


(A-36) 


Rewriting  Eq  (36)  with  a  dm  spherical  harmonic  subscript 
and  substituting  into  the  expanded  form  of  Eq  (A-30)  gives 


Weak-Form  of  the  Functional 


The  functional  whose  Euler  equation  is  the  even  parity 
Boltzmann  equation  of  Chapter  III  is 


where  the  inner  product  <f,g>  is  defined  as 


and  *  means  the  complex  conjugate. 

The  minimizing  function  of  this  functional  is  the 
function  I*  which  is  a  solution  to  the  second  order  even 
parity  Boltzmann  equation  (Ref  10:169)  Therefore,  a  solu¬ 
tion  to  the  even  parity  Boltzmann  equation  can  be  found 
by  minimizing  Eq  (B-1). 

Another  more  useful  formulation  of  this  problem  is  the 
weak  form.  This  weak  form  can  be  found  by  imposing  the 
condition  that  a  function  which  satisfies  the  natural  boundary 
conditions  Eqs  (41)  and  (42)  and  the  even  parity  Boltzmann 
euqation  must  also  be  a  minimum  of  the  functional  (B-1). 

Let  the  functional,  (Eq  (B-l),  have  a  minimum  at  f . 

Then,  for  ail  n  and  e  where  o  can  be  arbitrarily  small  and 
01  either  s ign 


(E-J^ 


where  ''i  and  n  are  real  functions  that  satisfy  the  boundary 
conditions . 


Expanding  Eq  (B-1)  in  H'  +  eq  gives 


(B-4) 

(B-5) 
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(A  7^^  /(XA-  y?Xiy  dX 

(B-7) 

(B-8) 

(B-9) 

(B-10) 

(B-11) 

-^((yt- vX^/r  J./- 

(B-12) 

(B-13) 

(B-1 

-se{<AX^^yd/^ 

(B-15) 
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(B-16) 
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noting  that  (B-4)  +  (B~8)  (B-12)  +  (B"14) 

+  (B-io'i  =  FC^) 

(B-19) 

6  b 


and  that  K®,K^,G®  and  are  self  adjoint  operators  (Ref 
10:174)j  where  if  L  is  self  adjoint  then 

=  <f(A\iscA)y 

-  (  B-2 0 ) 

*  means  the  complex  congugate.  Since  n  and  '1'  are  both 

real  functions  (B-5)  =  (B-6)  and  (B-9)  =  (B-10) .  Therefore, 


collecting  terms  and  simplifying 


with  both  k'^  and  G®  being  positive  definite  operators  then  the 

2  2  • 
e  term  of  (B-21)  must  also  be  positive.  Note  that  e  is  al¬ 
ways  positive.  Now  in  order  to  ensure  that  FCt  +  cri)  >  FCi) 
for  c  ^  0,  the  e  term  in  equation  (B-21)  must  be  positive  or 
zero.  But  e  can  be  of  either  sign  therefore  the  coefficient 
oi  c  must  be  zero,  that  is 


Eq  (B-22)  is  the  weak  or  Calerkin  form  of  the  second  order  even 
parity  Boltzmann  equation.  A  detailed  derivation  of  equation 
(B-22)  can  be  found  elsewhere  (Ret  5:57). 


Appendix  C 


Derivation  of  the  Weak  Form  From  the  Galerkin 
_ Method  of  Weighted  Residuals _ 


It  can  be  shown  that  solutions  to  the  second  order  forms 
of  the  Boltzmann  equation,  by  using  a  variational  principle  or 
the  method  of  weighted  residuals  are  equivalent,  A  proof  of 
this  equivalency  for  the  even  parity  equation  is  outlined  below. 
However,  a  similar  proof  can  also  be  extended  to  the  odd 
parity  second  order  equation. 

The  starting  point  of  this  proof  is  the  even  parity 
Boltzmann  equation 

and  vacuum  boundary  conditions 


(C-2) 

O 


-K(^i>)\Scr,rV  -  -i- 


(C-3) 


I  o  r  xL.  ■  Pu  y  O 


In  the  following  equations  the  v  and  9.  dependences  will  be 
omitted , 

If  a  trial  solution  'i  is  assumed  where  4'  is  a  linear 

combination  of  functions  such  that 
A' 


(C-4) 


then  the  Galerkin  method  of  weighted  residuals  requires 


Substituting  Eq  (C-9)  and  (C-10)  into  (C-6)  gives 

^-  d^iS^Lj^,?iyAh 
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term  A  of  Eq  (C-11)  can  be  rearranged  into 


(C-11) 
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Using  the  boundary  condition  of  Eq  (C-2)  and  (C-3)  ,  Eq  (C-12) 


can  be  written  as 
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Eq  (C-12)  can  also  be  written  as 


Y  \/d/?  J  /l<dAd  s 

dlfr 


(C-14) 


where  j,  moans  the  absolute  value.  Substituting  Eq  (C-14) 
into  (C-11)  gives 

/l(<n  -q.ddtvLyy  V, 
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Appendix  D 


Expansion  Properties  of  Spherical  Harmonics 


In  Chapters  III  and  IV  the  anpular  dependence  of  the 
trial  functions  and  cross  sections  was  expanded  in  spherical 
harmonics.  Because  of  the  two  dimensional  angular 
dependence  in  y  and  x>  the  requirement  that  the  expan¬ 

sion  functions  should  form  a  complete  set,  the  expansion 
is  presented  with  m  and  I  subscripts  as  follows  (Ref  6:608) 
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(D-1) 


z  z  ^ 
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where 


(D-3) 


y  (^,x)  - 


c 


(D-4) 
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(D-6) 


If  f(P,X)  is  even  in  the  an^le  X  then  the  expansion  must  also 
be  even  in  x,  and  therefore,  (D-4)  can  be  rewritten  as 


(D-7) 


where  the  odd  iSinx  term  is  omitted.  Also  from  Eq  (D~5) 


(D-8) 


Therefore  for  an  even  expansion  Eq  (D-3)  can  be  Virritten  as 
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For  uii  :ven  X  exnans  i  c' t  f(a,v)  Eq  (D-2)  can  therefore  be 


written  as 
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£‘=a  ^z-i 


rs  /^ 

where  P^(fi*i2^)  is  given  by  Eq  (A-24)  . 
Eq  (D-13)  can  be  expanded  to  give 


and  from  Eq  (D-5) 


,  ^  /rj-^  '  ,  ' 


Therefore  (D-14)  becomes 


0(A.A') 


.;.  A’zo  ffn^= 

whore  in  means  that  all  terms  w'ith  a  m  =  0  subscript  must 


bo  divided  b v  two. 

Since  the  scattering  cross-sec t ions  are  real  the  expan 
sion  of  Eq  (D-16)  must  also  be  real  and  with 
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Eq  (D-lo)  can  be  written  as 
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Note  that 


n(z-xO  'rri-:X.CxS7n'X^-r  Son.  (d-19) 


The  angles  y  and  x  are  shown  in  Fig.  2. 

Similar  derivations  would  give 
/  JL 


l--d7n.*^6  (n-21) 


Ry  inserting  Eq  (D-20)  into  Eq  (;\-22)  t'ne  even  parity  collision 
operator  can  now  be  written  as 


an.!  fr-'T.  Ecjs  (E-15)  and.  ('.>-17)  l!ie  invc'rse  odd  parity  collLson 
(.iperator  liq  (.A-^O)  ,  can  be  writlon  as 
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Appendix  E 

The  Synthesized  Boltzmann  Equation 

In  order  to  formulate  a  numerical  solution  to  the 
air-over-ground  problem  it  is  necessary  to  expand  the 
weak  form  of  the  even  parity  Boltzmann  equation  in  a  set  of 
trial  and  test  functions.  The  finite  element  space-angle 
synthesis  trial  and  weight  functions  of  Chapter  IV  will 
be  used  in  this  expansion.  Because  this  is  a  very 
tedious  and  lengthy  derivation,  the  more  obvious  albegraic 
steps  will  be  omitted.  The  starting  point  of  this  formula¬ 
tion  is  the  synthesized  even  parity  Boltzmann  equation,  Eq 
(58),  of  Chapter  III. 


and  they  are  given  by 
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Figure  9.  Surface  Normal  And  Particle 
Velocity  Direction  Vectors. 


and 


are  the  associated  Legendre  functions 


The  directional  derivative  is  given  as 

A 

A 


where 
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surfaces,  and 


I 


0 


i  o  r  the  vu  r  t  i  c a  i  (aid <? ■  a  f  t  !■■. e 
nr{''''1  era  cv  1  Lnds-r  ■  ' 


If  p., 


n  and  n  are  considered  to  be  unit  vectors  then 
r  z 


from  Fig.  7 


r  Bb^CF  ^1 

and  therefore  for  the  top  surface  of  the  cylinder 

and  for  the  bottom  surface 
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also  CD  =  HB  =  projection  of  p  unto  the  x-y  plane 
and  therefore, 

CD  =  HB  EB  SinO 


(E-12) 


=  -\/l-Cos^0  =  l-y^  (F-13) 

and 
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In  Appendix  D  the  even  and  odd  operators  were  given  as 


and 
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where  m  means  that  all  terms  with  a  m  =  0  index  must  be 
divided  by  two,  and 


Also  the  inner  product  is  defined  as 


(E-17) 


^ ^  forrealf  (E-18) 

Usins  Eqs  (E-16),  (E-15),  (E-14),  (E-11),  (E-12)  and  (E-6)  and 
noting  that  ,  Qn  ,  B-  and  B-  are  all  real  functions,  Eq  (E-1) 
can  be  expanded  to  give 
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where 
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and 


Eqs  (E-19)  to  (E-36)  is  an  expansion  of  the  first  term  of 
Eq  (E~l).  Eqs  (E~37)  and  (E-38)  is  an  expansion  of  the  second 
term.  Eq  (E-39)  is  an  expansion  of  the  third  terra.  Eqs 
(E“4U)  to  (E-46)  is  an  expansion  of  the  right  hand  side  of 
Eq  (E-1).  o^,  o  and  are  functions  of  z  and  they  must 

be  included  in  the  spatial  or  dv  integrals.  In  cylindrical 
geometry  with  azimuthal  symmetry 


A\/  -  2  TT/^cJ/'d'z 


(E-49) 


and  /ds  means  an  integration  over  the  surface  of  the  problem 
cylinder,  Fig.  9. 

for  the  air-over-ground  problem  with  an  exponentially 


■^3 


varying  air  density 


(E-50) 

(E-51) 

(E-52) 

where  0^(0),  a^(o) ,  o^^o)  are  cross-sections 

of  air  at  sea 

level . 

Z  =  the  height  above  sea-level 
sh  =  atmospheric  scale  height  ~  7km 
In  Eqs  (E-19)  to  (E-46)  the  integrals  are  separated 
in  the  space  and  angle  variables.  These  are  double  integrals 
in  space  and  angle.  However,  they  can  be  separated  into 
single  integrals  of  the  p,  x>  P  ^nd  z  variables. 


Appendix  F 


Angle  Integrals  of  the  Synthesized  Second 
_ Order  Boltzmann  Equation _ 


An  expansion  of  the  even  parity  Boltzmann  equation  has 
produced  twenty-eight  integral  terms  (Appendix  E) ,  By  a 
further  expansion  and  separation  of  the  integration  variables 
twenty  distinct  single  ancle  integrals  are  formed.  These 
angle  integrals  are  dependent  on  the  degree  of  the  spherical 
harmonic  trial  function  expansion  and  independent  of  the 
problem  parameters.  They  can  be  evaluated  once  and  thereafter 
used  as  a  part  of  the  problem  input  data.  In  this  research 
project  these  integrals  were  numerically  integrated  for  each 
combination  of  the  I,  m,  k  and  n  expansion  subscripts.  They 
were  then  stored  as  a  matrix,  and  selected  products  were  used 
to  produce  each  of  the  twenty-eight  angle  integrals  of  Eqs 
(E-19)  to  (E-46)  in  Appendix  E. 

These  twenty  integrals  are 
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A  numerical  evaluation  of  all  twenty  integrals  was 
carried  out  for  a  third  degree  spherical  harmonic  expansion. 
This  evaluation  showed  that  these  integrals  are  equal  to  zero 
for  many  combinations  of  the  Am  and  kn  subscripts. 

Integrals  (F-10)  to  (F-14)  are  a  part  of  the  surface  inte¬ 
gral  term  which  has  been  partitioned  into  the  outward 
(fi*h  >  0)  and  inward  (fi*n  <  0)  directions.  This  partitioning 
was  incorporated  into  the  weak  form  derivation  of  Appendix  C. 
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Appendix  G 


Space  Integrals  (Bicubic  Splines)  of  the  Synthesized 
Second  Order  Even  Parity  Boltzmann  Equation 


A  trial  function  expansion  of  the  spatial  flux  dependence, 
in  the  weak  form  of  the  even  parity  Boltzmann  equation,  has 
been  carried  out  in  Appendix  E.  Bicubic  polynomial  splines  in 
the  P  and  z  variables  were  used  to  form  a  tensor  product 
space.  These  splines  are  twice  continuously  differentiable  and 
have  non-zero  integrals  /g^B^(x)Bj  (x)dx  for  all  |i-j|>4. 

After  a  separation  of  the  p  and  z  variables  of  integration 
seventeen  distinct  integral  forms  are  produced.  These  integrals, 
which  include  the  source  integrals,  must  be  evaluated  over  the 
entire  problem  domain.  The  space  integrals  of  Appendix  E  are 
selected  products  of  the  following  seventeen  single  integrals. 
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The  Source  Integrals 
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cubic  polynomial  z-spline 
cubic  polynomial  p-spline 
atmospheric  scale  height 
outer  radius  of  the  problem  cylinder 
problem  cylinder  height 
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H(z)  and  H(p)  are  the  source  interpolating  functions  (linear 
Lagrange  polynomials). 
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An  Expansion  of  the  First  Scatter 
Source  in  Legendre  Polynomials 


In  Chapter  IV  the  first  scatter  source  was  defined  as 

-  <5Yz,A./i-)  (H-1) 

where 

({i^(r,z  ,^2  =  direct  fluence  of  Chapter  IV,  Eq  (74) 

a®(z,U*n')  =  scattering  cross-section 

The  usual  Legendre  polynomial  cross-section  expansion  will 
now  be  carried  out.  Also  the  even  and  odd  parity  first  scatter 
source  expressions  of  Chapter  IV  will  be  derived.  Expanding 
a®  in  Legtndre  polynomials  and  using  the  addition  theorem 
(see  Appendi'.  D) 
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where  m’  means  that  all  terms  with  a  m  =  0  subscript  must 
be  divided  by  two,  and 
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From  Fig.  5  and  Fig.  2  it  is  apparent  that  y  '  =  yd 


and  X  ^  =  0 
Therefore 
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and  using  the  identity  (Ref  18:96) 
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and  a  little  algebra,  the  even  and  odd  parity  first  scatter 
sources  can  be  written  as 
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A  Derivation  of  the  Total  Particle  Fluence 


In  Chapter  IV  a  trial  function  expansion  of  the  even 
parity  angular  particle  fluence  was  given  as 
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The  A-  .  „  mixing  coefficients  are  obtained  from  a  numerical 
1 ,  j  ,  K.  ,m  ® 

solution  of  the  second  order  synthesized  Boltzmann  equation  of 


Chapter  IV. 


The  angular  even  parity  fluence  is  also  defined  as 


An  integration  of  Eq  (1-4)  over  all  directions  gives  the 
total  even  parity  particle  fluence  4'(p,z),  and  also  the  total 
particle  fluence  tfo,;’!.  Ibis  is  because 


/iTT 


(1-5) 


Therefore 


^C^z)  -- 


Eq  <I-l)  will  now  be  integrated  to  give  the  total  particle 

auence  at  petition  (o.a).  Uaing  the  orthogonal  properties 
of  Legendre  polynomials  this  integration  is  carried 

follows. 

The  tero  order  associated  Legendre  Eonction  rs  defined 


as 


'PC^)  -  / 

'OjO 


(1-7) 


•  1  •  rn  ,'T-n  bv  Fq  (1-7)  and  integrating  over  all  o  and 

Multiolving  Fq  (It)  nq  i 
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where 


'2P' 


and  (see  Appendix  D) 


for  ra  5^  0 
^  for  ra  =  0 


(I-IO) 


-  2^  ^  ■=•  ^//^rT' 


(i-ll) 


therefore  Eq  (1-9)  becomes 


2  7r.z//^ 


(1-12) 


and  the  total  particle  fluence  is 


jrz  z/e 


where  2^  -^n  ~  O 

The  angular  particle  fluence  is  given  by 


A)  --  A)  -f-  X{/’Z,  A)  ( i-u) 


where  x(P>^>“0  is  the  odd  parity  fluence,  which,  is  defined 
in  terms  of  the  even  parity  fluence  and  source  by  the  follow¬ 
ing  expression 


( i  -  L  3  :i 


1 
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Appendix  J 
Computer  Subroutines 


A  computer  program  has  been  written  to  solve  the  synthe¬ 
sized  Boltzmann  equation  of  Chapter  IV.  This  program  is 
designed  to  use  a  trial  function  expansion  in  bicubic  splines 
and  spherical  harmonics,  and  to  perform  a  first  scatter  source 
interpolation  using  linear  Lagrange  polynomials.  The  program 
in  an  assemblage  of  several  subroutines  which  collectively 
perform  the  following  tasks. 

1.  Computes  all  single  space  and  ancle  integrals. 

2.  Combines  the  single  u  and  x  angle  integrals  for 
all  combinations  of  the  spherical  harmonic 
expansion  subscripts. 

3.  Combines  the  single  o  and  z  integrals  for  all 
coinbinat  ions  of  the  bicubic  spline  expansion 
subscripts. 

4.  Assembles  the  coefficient  problem  matrix. 

5.  Computes  and  interpolates  the  first  scatter 
source. 

6.  Assembles  the  source  vector. 

7.  Checks  for  symmetry  and  diagonal  dominance. 

8.  Solv'-'s  for  the  expansion  (mixing) 

coefficients  bv  the ’method  of  successive 
over-relaxation. 

9.  Solves  for  the  total  particle  fluence. 

A  ten  point  Kewton-Cotes  single  integration  routine  was 
useL;  t(.i  nume  r  i  ca  1  1 V  integral-'’  th  i  r  t\'-se  von  integrals 

of  foppendicos  F  and  G.  Idiis  integration  routin-c  is  an  in-'noi; 
suoroutine  of  t.ho  Air  Force  Ai.-ronan t  i ca  1  Svstenis  Hivision, 

'r  i  :h !  -  i’a  1 1  e  r  son  Air  Force  ‘'ase.  The  overall  program  logic 
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has  been  written  in  a  manner  whereby  this  integration  routine 
could  be  used.  The  program  is  written  in  Fortran  V. 


Listing  of  Problem  Subroutines 


PPDGPRM  MR  IN 

D I  MEN?  I  CL  0 :  3 «  U :  5;'  *  "  H  •  ‘ »  I  £>  »  DU 3  H  '  0 :  ?  -  0 : 

1  ■'  t  L  H  *  !*  ji' ^  I  *  L  '  Z'  ij  * 

rUMENCIDN  PNZ  ^  -  i :  7;.-  •  PflP  •  -£:?>  .  FMK  •  -E:  7>  .  IPZ  '  7^7. 

1 7 1>  7 1  R  :'  ?  XZ  E V L  t  0 :  4  >  «  ZVZ  •  7  ’  5  •  5>  '  SVX  *  7 «  5 «  3>  •  Z  VP  '  7  *  5 « 

1  O  ?  R  -  '.5 Oil  5  0)  !>  XSODL  *'0:  4 >  •  !  .’C  '-'0:  4  • 

CDMMDN  R 1  PN!X  !>  I  RE  R  .  H  •  ■  I  TYPE  •  I  T  7  •  I  TTY 
COMMnN.-B£  -  IZH  EE  -  IP 


-  PL ' E. E- S>  . DUPL  ' 

'  >  F'  P' 7  j  7 '  €• .  '  '  i-  F'  !  ' 
,  RPD  E)  -  EE  ■  E 


♦  THIS  POUT  I  ME  -OLVEE  THE  EVEN  PRPITY  RtaZOTPOPIC  EOLTZMRNN  ♦ 

♦  EC'URTION  RND  THE  RIP-OVEP  UPOUND  PPOBLEM  FOP  THE  PERU  T I  ON  ♦ 

♦  PRTE  FLUEriCE  BY  CRLLINU  R  NUMEEP  OF  EUBPOUTINEE  RND  UEINC  R  ♦ 

♦  TEN  POINT  NEUTOri  COTEE  INTEGPRTION  PDUTINE. 


P'ERD  '  ’  ♦ '  Er<D  =  1 E O  '  DZ '  ZE *  F'E  ?  DP' »  Z 1  ’  F'  1  *  R7  H  ?  L MR.’ 

PERD  '  ♦  *  END  =  1  E  i.i '  'r'D  ^  HI-  >  E  I  UT » I'l 
PPiriT-*.  PPOBLEM  INPUT  DRTR 
PPINT*. 

PPINT*.  DZ=  -DZ'-  DP=  -DP-"  PE=  «PE>  ZE='-ZE'  REH= 

1 »  RE  H  •  Z 1  =  ^  Z  1  »  P  1  =  ’  P  1  •  PELRI  ’RT  I  ON  FRl  TOP  =  '  '  I'l 

PPIt^T*, 

PPINT*.  -'  LEUENDPE  EXPRNEIDN  P  •  LMRX>  =  -  ■>  LMR;E 
PERD'  ♦’♦«END=1E0'  ':XEEVL'.I>  k  I  =  0»LMRX> 

PERD  '  END=1E''| .  -  XEDDL  '  I':'  '  I  =  0'LMRX> 

pcRii enD=  IE 0>  ■  ;EE  ’'I>  •  I  =  0-LMRX'> 

PPIr^T♦^  ■  YD=-  lYD'  HB=  •HB«  '  EIGT=  ’  .  EIGT 
PRINT*.-'  ■ 

PRINT*.  '  LEGENDRE  ODD  EXPRNEION  ECRTTEPING  CPOEE  EECTION.  RLL  ODD 
IMDNENTE  RPE  ZERO. 

-  -  ’  M  T  ♦  •  '  '  '  ~  E! ,  I  '  •  I  ^  L ' 

—  [  r! T  *  • 

FFINT*.  LEGENLFE  CI'D  PRPITE'  EE'PRN.IOr*  ECRTTEPING  C  PD  E  E  -  I  EC  T  I  Cl ' . 
ILL  'L'  EVEN  moment:  RPE  ZERO. 

PRINT*. 

PPINT*.  .  ■  .  XEDDL  •  I  '  .  I  =  ri,  LMR:  :  • 

PRINT*. 

PPINT*.  LEGENIPE  E,"PRNEIDri  rPCEE  EECTION  'ODD  RND  EVEN'. 

F'C  rriT,, 

-  I  r '  * 

i—  P-  ^  r  ^  •  * 

T  -  ■  L ^  i  ■  ♦  ■ ■  zi: 


•L  " 

rT-r  . 


J 


PRINT*.  THE  SPHEPICRL  hRPMDHIC  CDEFFICIENT-  'L.H'  , 

PRiriT*. 

no  5  L=u.LnR>: 

PRINT*.  '  .CL'L.N;:- .M=U.L:.' 

PRINT*.' 

PRINT*.- 
>■  1  =  ij . 

Xc:  =  c'*hTRN  ■:  1 .  :> 


TOL=. OuOOl 


f  DlJNT  =  E0  0 
ITVPE=1 

PRINT*.  RNbULRR  iriTEGRRLS  DF  IPHERICRL  HREMnNICS  FDR  THE 
PRINT*.  6RLERKIN  'IDLUTION  TD  THE  EVEN-PRR I T',-'  FDRN  DF  THE 
PRINT*.  ■  NEUTRDN  TRRNSPDRT  ERLIRTIDN. 


PRINT*.'  THEIE  RRE  INTEGPRLS  OVER  E*PHI  DF  THE  RERL  RNT  OR 
PRINT*."  PRPT  OF  EXP'-INX'.  EIRRNDEri  IN  SINE:  Rr^rl  COSINE:.'' 
PRINT*. 

no  15  ISH=1.1E 


L  HLL  S  F'HEF'  '■  LNR!'! .  i  H  '  U .  U .  I  i  H  '  .  !  ■  1 .  Xti .  TDL  ?  K  OLINT.' 
PRINT*.  '  HNHULRR'  'HRR'MONIC.'  INTEhR'RL='  .  I  "  H 


no  10  n=o.lnr:s 

F'  F'  I N  T  * .  '  '  .  -  H  '  M .  N .  I  H  - '  .  N  = .  L  M  H  ■ 

PRINT*.' 

PRINT*. 


CDNTINl.'E 
:a=-i .  0 
:  ■:s=  1 .  0 

ITYPE=£ 

PRINT*.  THE  RSSDCIRTEr  LENGENDRE  PDLVNONIRLS  INTEGRRLS  FOR 
PPIr^T*.  GhLEPl  IN  SOLUTION.  " 

CR'INT*.  THESE  RPE  THE  PDLVNDMIRLS  P'L'L.N"'  'RISDC.  LEGENDRE'' 
PPItJT*.  INTEGPRTED  OVER  THE  INTERVRL  '-1.+1‘. 

PRINT*. 

DO  £5  IP=1.P 

CRLL  POL'i' '  LMR::.  PL  •  1.  1.  IP>  jLI.LE.TOL.  FOUNT' 

PRINT*.  HSSOCIRTED  LEGENDRE  POL'i'NOMIRL  I NTEGRRL=  '  .  I P 
DO  £0  IM=1.KT 

!=-p'rf^T*.  ,  .  FL  IM.  i  T'  IF  '  •  IT- 1  . 1  T 

F':ir-T*- 

f '  P  I  '1  T  *  . 


CONTINUE 

CRLL  iNirLE  '  LMR!:.  lL*  .H*  PL.  LH.  .  IbT . !  !:,EVL.  X  -DDL.  FT  .  !■■!  _■  > 
ITT',  =1 

NF^'ZS-Sl'  DS-t-1.1 
Ii  0  :■  f  I  I  =  I ' .  r  F ' ,  + 

P'tJ'  '  '  T  '=  S  !■*■’!♦  FS 

-  '  ;  -  [  ■  r:’  ,  f  . 

' '  s :  = '  * 

[ TVPK= I 


EVEN 


THE 


rrr-.r*.  ..mI  I'Ct;  ~np  •'--i.Mi  r  ;  iTpr;t;'u.i_ 


>  I  r  ^  .  I 
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it:  =  i 

ppirtT*.  mhtpi::  vrlue:  fop  z-zprue  laiiPCE  integprlz 

ppIr^T♦. 

LD  EP  1=1 <3 

L  hLL  ;  DUPir  E  '  3  vz  -  nzs  -  tol  -  k  ouht  • 

CONTINUE 

n:  :=  -PE-pr'  •  DP+i .  1 
DO  6  0  I  =  C- N. 

PN' :  ■  I '  =p  1  + 1  ♦DP 

IF '  PN'; I ' .  EC.  0.  O  'Pn:  ■  •  i  .  =  i ,  oe-p 
PNP  ■  I  :■  =PNP  •:  I  ■ 

npp=n:  :+E 

NT  =  I  T  +  KT^  '  '  NZ  :-1  •  +NZI^-NP3-r:'  > 

PPINT^^ 

ppirfT^, 

PPINT^.  THE  TOTRL  NUNEEP  OF  TPIRL  FUNCTIDNZ  =  -NT 
PPINT^'  THE  NUNPEP  OF  LEGENDPE  FUIDITIDN:  =  .  1:  T 
PPINT^- 
FPINT^. 

PPINT^.  PPOELEM  REONETPV  INPUT  DRTR.  ■ 

PPINT^.  NDDE'Z'P'  ME3H  C ODPD I NRTE ' Z - P ■ 

PPINT^. 

DO  h5  mp=i;i, hZP-P 
DO  P*  M L  = '  I '  r p  _  —  !■ 

F'  F'  I  f  ^  T  ♦ «  '  >  M  F' «  «  'Ml*  '  '  *  P'  N  Z  '  f'1  F'  '  '  •  F'  N  P  '  T'l  C 

PPINT^. 

PPINT^- 

ITTPE=P 

ITI^O 

ppirp^^  NRTpr:  vrl'.e:  fop  p-pprce  itPECPRc:, 

FPir<T^, 

DO  pn  1  =  1,-; 

lrll  zflinE' : pp-np:- tol 'FOUNT • 

CCNTINIJE 

it:=e 

ppiNT^'  MRTPi::  vrlue:  fop  p-:prce  poupce  integprlp. 

PF INT*. 

T,  p  ^  1*1  —  1 

3  Hi  L  -Z',  F>  F  ■  ZvC  .  rP  :  ,  TOL  -  •  iPF-NT'- 

3  ZN  T  T  •  11  ip 

TO  COMPUTE  THE  I  l-TEF  PCLH  T  I  or  i  POINT  VRLUE: 

L  RL L  IN  F‘E L  T  -  Lt'lR! ' »  R*-  D ,  FT L.  >  F  f  iF' »  NP  _  '  NZ  .  '  .  I G  T »  HE'  V  D '  R  _  H  - 
TO  COMPUTE  THE  PP’OFLEM  ZOUPCE  '•■'EiZTOF 
RL  L.  EVE L  '  HF'  D '  -  *  .  VF' '  _  E  '  HZ  _  '  NF' '  LNR:  •  L  R  '  _  I  G T  *  i  T  <  N  T  ' 

TO  COMPUTE  THE  ZTIFFNECZ  RTPi::. 

>  hLL  .  T  I  FF  ' !_  h  '  F  F'  '  _  F'^  '  ri* _ '  UF'  _  '  I  gT  '  R  Z  '  FT  iF  '  F T^Z  '  L  MRL  '  N  T  - 

-r  'UJ-  f.-Hir-y 

i_  ■  'jp,-  V  ,  ,  y  X  . 

TO  ZOL'  E  THE  PPOELEM  MR‘"PI::. 

HZ.  L  '  OL  R  ',i .  ■  r  ,  M  ■  .  ■  • 

'0  I  'F  TPf-'  ; -L  v.  HLHw  Fli':.  HT  T'-JE  VCZrZ’. 

- l”^*. 

F  F  1  U  T  ♦  ,  1  Gp  T Z  Ti  -L  F  Hi,  I  HO  _  rt  /  t  E P  6 1.  -r -  f-  D  I  P  E  i  P  L J_3  ' . 

pcT-iT*,  -ruji'cr  cnTMT  "  .  "  ■  ~  i"  R  T  p  p  p  f.  r|ji' 

L  :•  I  '  -^L'  TO  T  -  ■■  r  I  P  T  ‘  -  1-iTT‘F  -  :-7  , 


loo 


MP= ' Mp;-? ■ ♦£ 

riz= '  zz-zi  '  Mz 
rip= '  pz-pi  '  rip 

DO  no  i  =  ow'iz 

DD  111]  J=0'MP 
p=pl+J^tip 

I F  '  P .  LT .  PMP  0  0 >  >  P=PHP  '  0:.> 
z=:Zi  +  nDZ 

L  HL L  TFLU'! '  F'NF' »  F'NZ ' !  !£■ »  F'  *  Z  *  f  T  ?  HZ  .  » f IP'  _  •  TF'H  I  «  DF'H  !  «  HB  >  .  I  b T  *  .  H '  'I'D  ' 

TFHI=DPHI+TFHI 

PPIHT*.' 

F'F'  I HT^ '  '  '  '  '  Z '  '  ^  F'  ^  »  TPH I  »  >  DF'H  I  • 

1  -TFHI 

CCr^TIHIJE 
GD  TO  no 

PPIHT*.  FPEHPTUPE  PPOELEM  BFiTh  END 

:tdf 

EriB 

:  UBF  OUT  I  HE  :  PHEP  ■  LMh:-:  .  BUSH .  K 1  -  KE  -  TOL  -  K  CUHT  :■ 

BiHErr  lor^  du;h  ■  o:  o:  s:> 

CCHHOH  B1  M.r<.  BE  I  SH  -  I  TYPE'  ITS  -  I  TTY 

♦  THi:  POlJTir<E  HUHEPICFiLLY  IHTEGPFtTEZ  THE  SIHES  FtHB  COSIHES  ZIHi'LE 

♦  Ir^TEGPHL;  OF  THE  RrH  ZOTPOPIC  EVEH  PRPITY  BOLTZHRHH  EC'ijHTIOn  B'  '  ♦ 

♦  UZE  OF  R  TEH  POIHT  HEHTOH-r GTES  IHTEGPRTIDr^  POUT  I  HE.  ♦ 


IF'  riH.LE.  S' GO  TO  EO 
IF  'IIH.EO.IO.'  THEH 
:  1  =  i:i . 

.'E  =  E.  ♦RTRH  '■  1  .  ' 

GO  TO  EO 

el:e  if  '  riH.EG.  ii;'  theh 
:  :e=e.  ♦rtrh  ■  i . 

: ;  1  ♦RT 'j'j '  ) ,  ' 

GO  TO  EO 

ELSE  IF  'ISH.EO. IE'  THEH 
: :  1  -t.^RTHH  '  1 .  ' 

,E-  -  '  !  . 

Zr^D  IF 

-pMTT'.i  ;p 

BO  m=o,lhr:-: 

iio  ji'i  lmr:' 

CRLL  mi^iRB  '  ."I  - :  :E'  TO!-'  BUSH  •  H  '  '  I  OUta ) 
I  F  '  B'J  :  H  '  H '  fi  '  .  L  T  ,  1  .  OE-  I  0  ■  BIJ  I H  -  H '  H  •  =  0  . 

LOrniMi  ,p 


'  !  pc  rir,  '  f  f  111-  "El  -  L  i  .  I  '  FE  •  r['^  ■  f  T'  F' 

r,  y  Mpfi  y  ;  Pf  i  '  I  'C  ■  -  ,  - 

::  •-  -•  p.i  ;  iT'.-F-nni' 


•  -u  '  r.  p' .  T  ‘  ,c;  T-it  --.'pppc-p  E-  -  p  JpTEr  LF'^E' 

♦  FOLvrnMiRL;  i.iHirn  “^ff  '^rft  np  the  s'.'Hthf:  izft  FV'^rt  ppptt' 


':Tr ;  n  ’  »J 
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1 


IF  '  IF.  LE.  6  '  bD  TO  ?ri 
IF  'IF.EQ.?:'  THEri 
:  :i  =  0. 

:  :e=  1 .  0 

GO  TO  50 

EL  IE  IF  'IP.  El?.  5'  TMEO 
' :  1  =  I? .  0 
:  :£=-! .  0 
EODIF 
CDNTir^lJE 
11,1=0 

DO  4  0  f=0,LMR:: 

DC  40  H=0»K 

IT=0 

Il,i=Il,i+l 

DO  4i‘i  l  =  o,lmh:: 

DO  40  M=o,L 
IT=IT+1 

CRLL  GUhD,:'!.'  TCL- D'JPL  ■  II.I-  it?  .1  DLifiT  ' 

I F  ■  D'JFL.  ,  1 1,1 ,  I  T  ,  .  L  T .  1  .  OE-  I  0  >  DUPL  •  I  i,l  -  I  T  •  =  0  . 

OOOTirtljE 
=  ETijpr< 

Er^^ 

pij'^oTicri  y,", 

[iiMEr'Tio.G 

rC'''POM  El  M.rt'EE  I:H  EE  IP  1:4  L  -  v  I  TYPE  •  I  T  ?  .  I  TTY 

ooMPOr,  Hi  pp:',i  h:h.  he  jp. jt  R4  i  ?.ih*lf 


*  THi:  FppcTior*  i:  urEL  by  the  teo  point  nemton  cote:  iNTEGERTior 

•  pCNTImE  'YijhT,  to  IELECT  the  INTEGFhTICN  FIjNCTIDN::.  it  ij:e:  the 

♦  IN'^rjPNHT  ION  IN  the  HEr'..E  CONNON  FLOCl  E  TO  DETEPM I NE  THE  FUNCTIE 

*  MHICH  I:  EEIf'G  INTEGFF-TED. 


IP  ■  ITT'.'.  E?  .  !  .GO  TO  5=^ 

IF'  IT',PE.EO..E''Gn  TO  5  0 

IP  ' I :h. eg. 1 ■  THEN 

','  =  r  n  *  ,  n  ■  ,  t'*' '  ■  *0  D  ' 


,  ,  ?  .  •  ■  0  •  .  C;  I  '  i'  ’*  '  '  -I ■*  I  I N  '  Y  '  *  ?  I N  '  !'!♦?,  '  '  ♦!,  0 1  '  fi*?. 

FETL'FN 

el:e  IP  '  I  :h.  E'' .  ?.  ^p .  I  :h. EC.  1 
','='101 ' '  ■ '  *1. 0 1 '  N-*'  '  c r '  H*-'  ’ 


0.  OF  .  I  :h.  EC.  1  1 .  OP.  I  ?H.  EC.  15  ■  THEfi 


V  =  CD;i  ♦CDS'  ♦  !  rr<  '  M^X:' 

RETURN 

ELSE  IF  ■:  ISH.  EQ.  S;>  THEN 

V  =  CDS  '.  rHX  :'  ♦"  ir<  ■  N^X':' 

RETURN 

ELSE  IF  'ISH.EU.9-'  THEN 

X  =  ■:  IS  □  S  :  S  >  ♦  IS  □  S'  N  ♦  X  '  -  N  ♦  S  I N  ■-  X  :■  ♦  S  I N  ■:  I'l  ♦  I'S ' '  ♦  -S:  I N  ■:  N  ♦ 

Erin  IF 

RETURN 

CONTINUE 

H=l-:  >♦£ 

IF'H.LE.  1.  ijE“1  I..I  '  H  =  U  . 

IF  'IP.EU.l':'  THEN 
V=H^PLF  'ST  -  K  .  rU  ♦PLF  >::X .  L  <  MS' 

RETURN 

ELSE  IF  ':'IP.EU.£>  THEN 

V='SC'PT  '' Fl '  ♦X^PLF  '  .X '  F  •  N  '  ♦PLF  '  .X»  L»  M-' 

RETURN 

ELSE  IF  IF.  EU.  S-:;'  THEN 

1'=  i.'RT  '  H  '  ♦F'LF  ‘X '  F  ?  N  '  ♦F'LF  '  X '  L  ?  M.-' 

RETURN 

ELSE  IF  'IP.EU.T'  THEN 

1'  =  !  '♦♦ci^F'LF  '  '.X  '  h  •  N  '  ♦F'LF  • ! !  •  L '  N  ' 

RETUFr^ 

ELSE  IF  ■  IF.  EU.  5.  CF.  IP,  EU.  r.  DR  .  IF.  EU.  S.'  '  THEfJ 
'1'=,' '  ♦F'L  F  '  '  rU'  ♦F'L  F  '  S’  •  L  «  M  ' 

RETUPri 

ELSE  IF  'IP.EU.t'  THEN 
V=PLF  '  S:.  f  .  N'i  ♦FLF  ■  SU  L*  N' 

ENr-'IF 

RETURN 

CCNTINUE 

IF' ITTRE.EU.E'UD  TC  SO 

IF '  ITS , EC .  1  ' UC  TC  !£0 

IF  ■  I  ,  EC  .  1  '  THEtf 

N=l 

N=1 

I  -  n  T  n 

I  E  !  “  I  .  S'"  .  S  'l-iRN 


ri=s 

UC  TD  rn 

ELSE  IF  'I.EU. S'  THEN 
C— s 
N  =  c: 

CD  TO  TO 


th 
vr.£,  rp- 
’  C'l  T  I  NT!^ 


I  n  i 


I .  S'' 


IF  ■  ITS.  EQ.  £;■  GD  TD  ISU 
IF  ‘I.EQ.l;'  THEM 
r'i=:-: 

M= 

GO  TO  7u 

ELSE  IF  'I.EO.c*  THEri 
N=:: 

b 0  TO  7  iJ 

ELSE  IF  'I. EG. S'  THEM 

M=1 

M=1 

b 0  TO  7  iJ 

ELSE  IF  ■I,EQ.4>  THEM 

M  =  3 

M=1 

GO  TO  70 

ELSE  IF  '^I.EiS.^;'  THEM 

M=1 

M=1 

b 0  TO  7 1 1 

ELSE  IF  I.EQ.G>  THEM 

M=1 

M=1 

Er^DIF 

:Pl  =  SFF'^^^■  .pm::'  jf-S'  .pm:: '  jp-e':.  ,pm:7':jp-i  '  .pm::'  jP' 

SF''i'='F'F  '  fM  f  —  JT  .  F'M:  :  '  JF'+ JT— .  FTL:: 'i  JF'+ JT~c:,'  .  F'M:  :  '  JF'+ JT— 1  '  .  F’M:7 '■  JF'+JT 

1 ' 

IF' IT'/PE.EC'.E'GD  to  100 
IF  'I.LE.S'  THErt 
V=E.'  'P  ' : .  H S  H  '  ♦  S  F'  1  ♦ :  F'E 
PETUFH 

ELSE  IP  'I.  EC'.  4'  THEM 
V  =  E:  'P  ■  :  HSH  '  ♦SPl^SFE 

PETIJP7( 

ELSE  IF  'I.  EC'.':.'  THEr^ 

S  P  1  ♦  ■  PE 

;  ITT'  .c  ri 
■scriTirr^E 

IF  'I.LE.S'  THErt 

IP  ' : Eo.  n.  •::=!.  me-S''i 

■i  =  :fi*:pe 

pETHPr. 

ELSE  IP  '  r .  Eb.  4.  OF  .  !  .  EG.  S  '  THEft 


SC  S  S'  I  .  S  '  .  -  ■  S-pEf' 
-  r  ;  *  ■  t  -  • 

.  =  ~ .  : c  r; 

-  ‘  ,  r.  r  z 

T  F  .  r  ,  M 
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THEN 


H=£ 

60  TO  140 
ELSE  IF  a.Eij.3> 
r^=i 

60  TO  140 
Er^rlIF 
CONTINUE 
IF  a.EC'.i:>  THEN 
N=? 

60  TO  140 
ELCE  IF  a.EO.cl.DP.  I.EQ.  5>  THEN 
N=1 
END  IF 

S  PH  1  =  :  PF  ■  N -  K  C -  PNX  ■:LP-3:>  •  PNX  -  LP-S^  .  PNX  -  LP-  1  i.  PNX  ■  LP:'  -  J-T' 

-F'Hti^HIF  '■  IH'  F'NX  '-.L.'  ^  PNX  '  L  + 1  '  '  .'-vJ 

IF'ITS.Eu.c'/'EO  TO  16  0 

IF  '  I .  EC.  1 .  OP.  I .  EC.  £;>  THEN 

'r' = S  F'  H  1  ♦  :  F'  H  3  ♦  E :  ■ '  P  >^1  ■-  H  3  H  > 

PETIJFN 

EL  IE  IF  'I. EC.  3'  THEN 

V=IPH1#:FH3 

PE TUP N 

EriDIF 

continue 

IF  ■:  I .  EC.  1  .  OF  .  I  .  EC.  E;:'  THEN 

'(■'=:  PH  i*:ph3 

PETU'^'N 

EL  IE  IF  a.  EC.?'  THEN 
'i'=:ph1'»?fh3*::. 

END  IF 

PETUPN 

END 

FUNCTION  PLF  '.X-  I «  J > 


♦  THi:  POUTINE  IELECT:  hND  COMPUTE?  THE  H??DCIhTED  LE6ENDPE  ♦ 

♦  FUr^CTION  liiHICM  I?  FEIN6  INTE6PHTED. 


.E.  .  .  '  E-!  i:' 


»  I  , 


IF  ' I . EC . 0 . HND . J . EC . 0 ■  THEN 

PLF=^1  .  0 

PETUFN 

ELIE  IF  '  I .  EC.  1 .  HND.  EC.  0  '  THEN 
PLF^: 

C'CTCcr' 

?L  ?  IF  I  .  C '  .  t  .  HN D I  ■  ’’’HFl * 


c  p  Ti  it  ri 


:?  I""  I .  .  .1.  -r'D. 

-  i  *  •  •  -  +  a  *  f;  ■ 


ELIE  IF  '  I  .  E .  c  .  HN  !■ .  J .  1 1*' .  1  '  THEN 


.  .-H'D. 


loo 


PLF  =  ::'^H 
PE TURN 

EL  -  E  IF  '  I .  EiT' .  .  FtNH .  J .  Eu .  lY.  THEf^ 

PLF  =  '  4c‘ .  ♦  ' !  .  ♦X^B  '■  ■  * 4b . 

PETIJPH 

EL-IE  IF  -  I .  EO.  ?.  FND.  J.  EC.  i:>  THEN 
PLF  =  ~  '  Z  C'PT  '  H  '  ♦  '  iibb'^X ♦♦!=.■+ 7c.' 4b 
PET  UP  M 

EL'E  IF  ■  I.EO.X.RtiD.  J.Eri.c.  THEN 

PLF=15.  ♦X^Fl 

RETURN 

ELIE  IF  '  I.EO.  b.HNIi.  J.EO.  THEU 
F'LF=-15.  ♦  'H**!  .  5  • 

END  IF 

PETUPN 

EUD 

rUBPDUTIHE  hHGLE  •  LMh>;  -  C  L  -  7 H .  PL  -  C  H  .  '  IbT. : !  lEVL  -  X  IDDL*  h  T  -  X  '  ■ 

B  I  r'Ert:  I CU  C  L  •  U :  ;  <  O :  X  ■  «  I  H  ■  O ^  U :  3 ,  1  3  .  .  PL  '  x  •  3  ’  3:  '  •  il  H  ■  3  -  3  -  3  L''  '  -  !  X  E'.'L 
14*'  ’ .  DLL  *  U  J  4  *  '  /.  _  *  H  t  4  * 


♦  THi:  PGUTIfC  R'IERELE:  the  TDThL  hH'3LE  IHTEGIX'hLI  fop  hll 

♦  coNBiuHTiorv:  the  iphepicFiL  hhpmdnic  tpihl  hhp  i.ieight  FurMiT: 

♦  iubicpipt:. 


PPirJT*.  TOThL  ir'TEGPRL  VRLUE:  iri  MHTPIX  FCPM.  OF  THE  IPHEP IC*^-. 

C'Pir;T»,  HHcr.trr.n:  tpihl  hub  '.iEIGHT  f^i.bJrTICU:  iJIEB  ir^  THE 

ppTfjT*,  GPLEP»  IH  rOLUTIOfi  OF  THE  EVEf^-HPC  I T','  BDLTZMHrir<  EOUHTIC’* 

PPIHT*,  the:e  vmlue:  hpe  ‘electeb  ppdti.'lt'  df  the  hizocipteb 
PPIHT*.  LEL'Er^BPE  FOLYHOMIHl:  hub  Hrjpp.^pp  hjtegfhl;  UHICH  i.EPE 
PPIHT*.  rOMPUTEB  EhPLIEP.  THF','  PEPPElEfiT  THE  rCHTTEPIf'G  hHB 
PPjht*,  r^^r^-:^HTTFF  Ir^G  hHGULHP  IHTEGPHI..  ^EFtJEL:  i.*HIlH  hpe  FCPr*'- 

PPJ^JT*.  hheh  :p'HEpichl  HHEMcriic:  hfe  u:eb  h:  the  hhih.ilhf  tpih^ 

PFIHT*-  hHB  i.*EIGHT  FUrtCTICH"  IH  THE  GHLEFlIf^  rOLUTICU. 
pp  : 

ppiriT*. 

BO  c'c’O  1  =  1 'E'X 

Il,i=l  . 

IP  •  P'" .  ;  thfh 


GO  re  1 1.* 

ELIE  IF  'I.EO.c'*  THEt< 
I  :h=3 
IPh=  1 
GO  TO  1  I i 

ELIE  IF  '  I.P'.  '  TMC-. 


-■’0  TO  1  *■• 
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.50  TO  1  0 

EL'.e  IF  'I.EC'.E'  THEN 

IPh=£ 

GO  TO  10 

ELSE  IF  kI.E0.7-<  THEN 
ISFi  =  S: 

IPh  =  E 
GO  TO  1 0 

ELSE  IF  <  I.  EG.  S'.  THEN 

ISH=5 

IPh==' 

GO  TO  ?E 

ELSE  IF  .I.EG.G'.  THEN 
ISH=P 
I  PH  =  4 

GO  TO  1 0 

ELSE  IF  'I.EG.l?.  THEN 
ISh  =  E 
IPH  =  -r. 

GO  TO  1 u 
END  IF 
GO  TO  SG 
DO  SO  i=o«lmh:. 

DO  so  N=0«I 
IT  =  ri 
11,1=11,1+1 

DO  3  I'l  L=0'LMP! 

PD  SO  ■■1=0.  L 
IT=IT+i 

L  H  '  I  '.I  •  I  "1  •  I  '  =S  L  '  1  .  tf ,  ♦r  L  '  L  ^  M  '  *  3  H  ,  M .  ri .  I  .  H .,  ♦F  L  ,  I M  •  IT.  I F  H  ' 
CONTIN'.iE 
GO  TO  SS 
DO  SS  f  =n.LNh: 

DO  :S  N='>.f 
IT=0 
iM=ri,.+  i 


:t-=  :  1 

H  '  1 1"  .  1  T  ’ 
OONTlrilJE 
IF .  I  .LE.G 
IF  'I.Fr,. 
DO  4i‘  •  =0 

jiH  4ii  N=r, 


SO  4I.  S-l.: 


H  ■  f  ^  ^ '  i  . .  H  '  ♦  ^  L. 


[T-  ■ 


OP.  I  .EG.  l'--GC  TO  cl'"-. 
:l.  OF.  I.  EG.  S'-'  THEIi 
LNi- 


■PMF  t  =  .  H  ■  r’.  '  ,  i  ^  •  .  ■  .  ■ 

■  M  '  I  "  I  '  I  '  ='*1  '  ^  .  H  ,  *r !  .  , 


*  ,  Tirt'C  1  *p|  ,  II,;  ,  I  T  , 


'.n  r'~  -II' 


[ .  S','’. 


1U7 


TEMP  1  =PL  '  I M  -IT,  7::'  +PL  a  l.l  -  I T  -  3> 

ch  '  I  u .  I T .  I  j  =i:l  (k  !.  ro  ♦cl  < l  •  m  -  ♦temp  i  ♦  i  h  '.:m  .  m  .  C  ' 

EMDIF 
CONTINUE 
GO  TO  £00 

ELCE  IF  '/I.EQ.LE::'  THEN 
IxR  =  3 
IPh  =  3 
GO  TO  45 

ELCE  IF  a.  EC.  £5:'  THEN 

I  IR=S 

IPR=£ 

GO  TO  45 

ELCE  IF  ■I.EC.£4.  THEN 
I  CR=t. 

IPR=5 

GO  TO  45 

END  IF 

GO  TO  55 

DO  5  0  h  =0'LMR:: 

DO  5  0  N=0<f 

IT=0 

IM=IM+1 

DO  5  0  L  =  ri>LMR!' 

DO  5n  M=0-L 
IT^IT+1 

m:=l-m 

R= ' 1 . - ' -1 . ' ♦♦m: ' 

I F  '  f't .  E I. ' .  i.i '  H  =  M  ^ . 

C  R  '  I  I'l '  I  T »  I  j  =r L  ‘  L  *  M  '  ♦♦£♦•'0  '  T  "  N  '  ♦:£  H  'T'l  <  TH  I R.'  ♦F  L  '  I  FI »  IT?  I  PR  '  ♦R^ 

IF'  I.GE.££.HND.  I.LE.£4:'G0  TO  £00 

CONTINUE 

IF  '[.EC.ER.'  THEff 
DO  c  0  f  =0-  lmr:  : 

DO  Ri'i  N=0«K 
IT=0 
T  1,1=11,1+1 
'  C  -  ij 

C  G  ^  L. 

I  T  =  I  ^  +  1 

m:==l-m 

R= '  1 .  +  ' -1 .  ' ♦♦Ml  ■ 

IF ' M. EC. 0 ' R=R  £. 

CR  '  I  FI.  IT.  I  :■  =CL  '  L.  M  ■  ♦♦C  +  CL  -  I  .  N  ■  ♦  i  H  '  M  .  N .  +•  '  ♦F  L  '  1 1 1 .  IT.-:.'  ♦R^\  I  ■  L 
GO  TO  CFii.- 
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THEN 


ISri  =  T’ 

IPFI=3 
IPB=3 

hq  to  1 

EL'-E  if  ’I.EO.lc:' 

I3R=I: 

I-B  =  ':- 
I  iri  =  F 
i:e=s 

IPH=':' 

IPB=3 
GD  TO  13U 

else  IF  'I.EU.l^-  THEff 
I;::h=5 
ISB  =  3 
IiB-3 
ISE=7 
IPH=3 
IPB=3 
GO  TO  130 
EL'E  IF  'I 
I  :h=5 

i:b=5 

I3D  =  - 

i:e='?' 

iPF=:-' 

IPB  =  T 
GO  TO  170 

el:e  if  'I 

I  :P  =  ‘:' 


EQ.14.  THEtl 


EO.  ll' 


THEri 


I  :b=- 
I  rr'^- 
i'E^'- 
IPH-'"' 
TPB-- 


i:b=I' 

i:b=3 
I  :e=7 
TPH=: 
IFB=3 
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I  :E=t 
r::ri=:i; 

IFH^5 
IPE=5 
GD  TO  130 

EL-E  IF  a.EO.EO:'  THEM 
I  Ih  =  F 
I3E  =  F. 

I  :ri=3 
HE =3 
ifh=f 

IFE  =  F 

GO  TO  130 

EODIF 

GO  TO  150 

DO  145  i=o-lmh::: 

DO  145  N=0-1 
IT=0 
11,1=11,1+1 

EO  145  L^O'LMh;- 
EO  145  M=0'L 
IF=0 
IT=IT+1 

L  H  '  I  M  ^  I  T  •  I  '  —  . 

EO  145  jF  =  0HriH:-: 

TEMP=0. 

:t  ij=:  ''OIL  '  jF  '  ■ '  igt-: toel  ■  jf  •  • 

IF  '  I  .  EG.  EO;. ! 4  lj  =  : E VL  ■  JF  ■ 

ED  140  j:=o<jp 

IF-IF+1 

Ih-E 

IF'  .1 1  •  E '  j  .  Cl I F  =  I F  ■  E 

TEMP'=IF*i.L'l  'I’i'^OL'L*  to  ■  ♦  i.  L  '  .i  F  «  J  .  '  ♦  ♦  c  ♦  '  .  H  '  _i .  '  f  i »  I  .  F  '  ♦  .  H  '  j  _  <  f'1 '  I  .  E'  ' 

1 j:  .  r'.  I  :e  ■  ♦:h  '  j:  .  n- 1  ie;.  ■  ♦fl  -  it-  if.  iff  ■  ♦fl  .  im.  ip.  ifej  +temp 
rorHiouE 

CF  .  I'.l.  IT,  I  .  =TEMP*:  4  IJ  +  GF  '  II, '.  IT*  I  ' 


I  :f=i 
I  :e=+. 
I  :E  =  r 


I  :e  =  ': 


ri_'p  re  ^  ■'■i-CM 


I  :e=  ■ 


el:e  if  •i.eo.eF'  theh 
I  !h=E 
I  :i'=E 
I  :ri=5 
I  :e=e 

IPR^S 
IFF  =6 
EriDIF 

DO  17S  r  =0,  lrr:-: 

DO  17'5  r^=o^^ 

iT=i;i 

IM=IM+1 

DD  17'?  L=0-Lr'IR': 

DD  177  H=iuL 

ip=i:i 

IT=IT+1 


L  H  '  I  liU  17’  I  '  —  U . 

Dn  17'?.  jf-M’LMr:; 

TEMF-m. 

::i  u=:::gdL'  jp  •  ■ :  irt-: ::ar'L  ■  jp  •  ■ 

DC  170  j:=o<jp 
rp=iF+i 

R=E* ' 1 . - ' -1 . ' •♦R? ' 

I F ' J : . ER . 0 . DP . M . EO , 0 ■ R=R  £ 

terc^i'l  '  L<  M '  *rL '  y . H .  *\:i_ .  jp,  j:  :h  ■  ji .  h«  ior  ■  ♦ih  >  ji .  ri  - 1  :e 

I '  fR  I  :d  ’  *:r ■  o: •  r-R  r  :e ■  •  ■*fl  -  it  -  if*  ipr  -  ♦pl  •  im.  ip.  ip?  '  ♦r+temp 
ccr^Ti'-ijE 

rp.  IM.  IT.  I  ■=rR.  Il,i,  IT.  I,'+TEMP*71  U*K?  -L.' 

LDRTir^ljE 

FPiru*.  MRTPi::  vrluez  fop  rurle  irtei3prl=  .  i 

DO  £l  0  I(,i=  1  .  f  T 

ppira*.  '  .  rp .  11,1,  IT.  I  ■ .  iT=i . f  T  ' 

PPIRT*. 

PPIRT*. 

l:□r^TIr^uE 

PETI.'PM 


♦  THi:  PCIJTIRE  compute:  the  IPHEPICRL  HRPr''OHlL  E.'.'PRHiIOr' 

♦  lCeffic  ieht:  op  r-oPHRL iertioh  roNiTpr.T-. 

PERL  IFF. IFF 


4 


'-F  Tn 


IDF^IDF^I 
FT=ir<F  IDF 
Gn  TO  100 
FT=1  . 

CL  '  LW'1  '  =:0PT  ■  ■  L^L+l  ■  ■  4^  ■  4*hT>=.n  .  j  .  ,  ,  ,  t 

c□r^TIr^lJE 

PETUPH 

EhD 

suFFDUTrr-iE  : 'PLiNE  ■  :p.”. h::: .  tol. f  dldt  ■ 

D I  MErt:  I DH  Fh' :  •  -e  :  7  •  • :  p:  :  ■  r .  7 .  c  ■ 

CnMMDH  -  Hi  ■  PHIG  I  Hi:  JF  -  JT  E4  L  » 1 


♦  THIi  POUTIHE  LDHPUTEi  THE  PICUI^IC  IFLIt^E  IHTEhF'hLI  E'i'  CRLLIHb  ♦♦ 

♦  H  TEH  POINT  HEMTCH-COTE:  IHTEGPhTIDH  POUTIHE. 


DO  '?  U  .1 F'  —  1  » r ^ . 

DO  GO  jr  =  i.ri:::i 
1:  F' ! ! '  J  F'  ^  J  i_  '  I . '  =  I ' .  ij 
JT  =  JC-JP 

I  F  I  H '  ._l T  ‘  T  .  I-  '  H D  TO  G 

I  F  =  4 

IF  LJT.EG.O'  THEli 

II  =  1 

GO  TO  5 

EL:E  if  'JT.EO.l'  THEH 
1  l=i 
GO  TO  j 

EL:E  if  '.'T.EO.i'  THEH 
^  1  =  ' 

GO  TO  ? 

el:e  if  'JT.eg.:;;'  theh 

1  1=4 

GO  TO 
EHDIF 
^  1  =  1 

IF  'JT.EO.-I'  THEH 
1  F=  ? 


GO  TO  5 

EL.E  IF  'JT.EO.-:;.  THEH 

LF=1 

GO  TO  S 

EliDIF 


r  or<  r  i  r4  :e 


1 

i  ^ 


1'  1  =- 

IF  '.IF.E0.4'  THEN 
IF  'IP.EO.l'  THEH 

;ff=z 

RETUPn 

el:e  if  'If.eo.E'  theh 
;ff=zi 

FETUFff 

EL  IE  IF  'IF.EO.?’  THEr^ 

:ff=z+:'>zi 

FETL’Ff^ 

Er<riiF 

EMI'  I F 

'i'  =  Z-J*E 

V  1  =  t  E  .  ♦E:*-*E  +  Z1 
IF  'JP.EO.::'  ThER 
IF  'IP.EO.l.'  THEM 

:ff=',' 

FETi.iPM 

EL:E  IF  'IF.EO.E,'  THEM 

:ff=vi 

PETijpn 

el:e  if  'If.fm.?.  thf.m 

-P'P=',  +  1 


£m:  :f 

'.,.'  =  V  +  F  ♦!'  *♦  - 
!  =v  1  -  !  .  ♦(!  ♦  *E' 

IF  '  pn.  E  ■  THErJ 
I!^  ■[p.r:0.l'  THEM 


err 
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D=::i-: : 
l,l=V-4^D^*  5 
Ml  =V1  +  1£. 

IF  ■JP.EO.i:-  THEIJ 
IF  ‘IP.EO.l.;'  THEM 

:pf=i.i 

PETUPM 

el:e  if  'Ip.eo.E'  them 

;pF=i.ii 

PETUFM 

EL:E  if  'IP.EO.:::-  THErJ 

;  P  P  =  I,1  +  :  >1,11 

Er<DIF 

EMDIF 

PETUPM 

EMIi 

lUBPDUTiriE  IDUPCE'  I V: !  -  M!  M  .  TOL  .  1  DUMT  ■ 
r  I  MEM :  I  cr^  :  v:-:  ■  r .  5  ^  i  ■ .  pm:  :  •  -e  ••  r  • 
roMMOru  Hi  Fr<; I  h4>:,iH'LP  fj  L-f 


♦  thi:  poutime  :elect:  hmd  compute;  the  iMTEFPnLHTEr  idupce  * 

♦  If^TEl3PHL  E'i'  C>^L.Lir<H  H  TEM  FOIr^T  MEI.iTDM  CDTEI  IMTEGFHTIOM  PDL'Tir-E 


Fn  Pm  lp=i.m::: 

DG  PM  LC  =  i-rf'  :-E 


'  L  E'  ^  EL  '  I ~  «  1* 

LV=LC-LF 

IF'LM.GT.  l.DP.L'.'.LT.-I:'GO  TO  PO 
f  F  =  4 

IF  'L'l'.EC.l'  THEri 
^  1=4 

GO  TO  1 M 


el:e  if 'Lm.E'm.  m.  them 

^  I  =  :: 


GG  TG  1  M 
Er^PiP 

if  'L'i.Eu.-1'  THEI! 


EL;E  if  'LV.E'M.-c'  THtti 

f  1  =  1 

f  F=^ 

C-C  TG  1  M 

EL  _  E  IF  '  L'l  .  E'"' .  -  3  '  THET' 


TU=  1 

;  r  F 

-  .  r;  _  TT  ^ 

L  E'  “  L  f’  ~  t  * !  ,  ^  _ 

.-.I  I  ■  M  ■  ,  .  '  r*  f 

U4 


SVX '  LP.  LC.  I  -' =;VM'LP*LC  .  I  '+TV 
IH=IH+1 
c  OHT I  r<i_;E 

PPIMT^.'  NhTPIX  SCUPCE  IHTEhPRL^'.I 
DO  9  0  NP=1^N/::: 

F'P  I '  '  '  ?  1 V’ !  ‘  MP »  MC  !■  I ?  Ml  =  1  >  WKZ  ~c!.‘ 

PPINT*. 

PETUPM 

EMD 

Fur^cTI□r^  HIP'  iH.xi.x£.:;> 


♦  THIL  RDUTIME  COMPUTES  THE  LIMERP  LRGPRMGE  POLYMOMIhLS  UHICH  ♦♦♦ 

♦  RPE  EEINb  IMTEGPRTED  RS  R  PRPT  OF  THE  SOUPCE  IMTEGPRLI. 

IF  'IH.EO.S,:'  THEri 
HiF= ' 

PETUPM 

ELSE  IF  'IH.EU.l:'  THEM 
HIF=  'X-Xl)  vXE-:  'i:' 

EMDIF 

PETUPM 

Er^n 

j  UEPOUT  I  TIE  D  I E  El  T  '  LMR!  M  RFTi  <  F  MZ '  PMP « tIF'  I  •  MZ  I  '  ”  I  G  T  •  HE  »  VD  *  R  I H  ' 

D I  MEt‘:  I OM  RPP  ■  5  •  5 .  J  ■  •  PMZ  '  -Z :  ?  ■  •  PMF  -Z :  7  • 


♦  TMi:  PCUTIfiE  CCMPUTES  THE  DIPECT  SOUPCE  POIMT  VRLUE:  UHICH  RPE  • 

♦  UCEr  IM  THE  FIPTT  CDLLI-IOri  CDUPCE  I MTEPPCLRT I  OM, 


PPIMT*.' 

PPIr^T♦.  IMTEPPOLRTIOM  POIMT  VRLUE L  .  ELEMEMT  ■:  I  -  J> 

1  ,P'i 

PPIMT^.  ■ 

I  (.1=0 

Eo  SO  l=o«lmr:-: 

DO  so  M=0-L 
11,1=11,1+1 

DC -HO  IT=i.MZS-Z 

:c  -io  i^iWiFz-s 
ZH--HLZ  '  I  i  ' 

P  T'rME',  I  —  i  ‘  •♦cl  +  _H  +  +c  • 

I  F  '  F' ^ H  .  E L  ,  I J  '  I ^ 0  TO  G  U 
UD=ZH  pzh 

IF ' UD. EO. 0 ' GO  TO  ZO 

T  H  l_i  =  _  I  L""  T  ♦  '  E ' !  F'  '  —  H  £'  H  _  H  ’  —  E ! !  F'  '  —  F'  f '  Z  '  I  T  —  1  '  H  .  H  '  '  •  R  .  H  LI  D 
GO  TO  ;0 

+RI  ■-  '  1 9  ~  •F'  'c  ,  -c-ri-  ,  T  T_  <  .  u j '  u  ■  ■»£  ric- .  T  -  t  . 

1  ^  •  J  •  Z  '  -  •  T'  * 0.'  F' '  —  ^''•i  ■  •  .  ^ r •-‘-I'l  ■  '  .  * F  ^ - 1  H 

1  r  L_  '  [  T  - 1  ■  -  H  ’ 

’’ZMT  [‘"Ip 

-t "  T  M  T  * ,  '  M  r  9  1 1  ’  L  L’  E  '  F  i":  -■  I M  r  ~  -  F  rj 7  '  Z  ‘  ^  S  I  ^  z  '  ^  U 


-  F  [  r  ■  r  ♦ . 

[iP  a'Z’  r'P  =  i.PZ'-i. 

-  C  '  r  ►  1  r  *  .  ,  .  r- 


FOP  POIMT 


r-'U  .  . 


b 0  TO  7  i.i 

PRINTS.  IMPROPER  IMPUT  DhTh.H  NODE  CRHNOT  EE  hT  THE  EUPIT  POir-T-Z 
1=  .PNZ'IT-l' 

7  TOR¬ 
RE  TURN 
END 

VUE  ROUT  I NE  EUEC  ■  FiPD  -  VVZ  ^  'i  VP  -  VE  -  NZD  -  NR  7  ,  LNRK  -  C  R  -  7  I GT  .  T  T  •  NT'' 

D  I M  E  f  I  _  ION  R  f"'  D  '  7  ’  ’  V- '  ^  V  V  Z  •  7  ^  7  *  ^  V  V  R'  '  7 »  5  !<  V .  '  ’  E-  ■  7*  iJ  '  '  L  R  •  1  V  ’  ■'  iJ  ' 


♦  THI7  ROUTINE  R77EMELE7  THE  PROELEM  ZOUPCE  VECTOR. 


PR  I  NT*-  ■ 

PRINT*-  ZOURCE  MRTRIX  CDLUMN  VECTOR.- 
PRINT*- ■  ■ 

17=0 

DO  90  JZ=1-NZZ 
D  0  '?  IJ  -I  R-  =  1  -  R'  V 
JR=U 

DO  :3  0  jl=o,lnr:: 

DO  SO  JN  =  I) -JL 
JR=  JR+ 1 
I  Z  =  I  Z  J-l 
7  E  ■  I  Z  '  =  I' . 

DO  77  N=l-7 
Z  C  =  1  ,  I J 
TEMP  7  =  0. 

I R  =  iJ 

DO  7  0  IL=0.LMR'': 

DO  70  IN=0- IL 
IR= I R+ 1 
TEMPS- 0. 

DO  CO  iz=i-rcz-s 
DO  CO  ip=i-nr:-s 

n:  P=IR-JR 

IP' rCZ.CT, l.OR. irZ.LT.-S'CD  TO  CO 
IP'  ICP.CT.  l.GR.  ICP.LT.-Z. 'CD  TO  CO 


fi‘^=  17 
GO  TO  '^1' 

EL:E  11^  'N.RU.Z-  THEN 
NZ=  1 
■NR  =  1 

Z  ~  ^  Lj  C  ^  . 


LU) 


r  — 


'  rv  F'"'.  4  . 


TL.cri 


150  TO  5  0 

ELSE  IF  (H.EC'.S::'  THEH 

hZ=  1 

hP=E' 

NR  =  £6 
bO  TO  50 

ELSE  IF  ':N.EQ.6)  THEN 
NZ=E' 

NP  =  3 
r^H=E•7 
GO  TO  50 

ELSE  IF  ai.EiS.?:'  THEN 
NZ=S: 

NP=:; 

NR  =  £3 
END  IF 

TENF'  1  =  S  VZ  '  JZ ’  IZ^NZ  '-*-  '.•'P  '  JF  «  IP»  ffP >  ♦RPD  I Z •  IP  •  I R  '  ♦)_  H  '  JR »  I R '  NR.' 
TEMP£  =  TEMPl-t-TENP£ 

COraiNiJE 

TEnP3  =  TEMP£  +  TENP:: 

IF  '  NR.  EC.  £3.  DP.  NR.  EC.  ^6'•  Si:=-1 .  0 
IF 'N. EC. 7) GO  TO  7£ 

TENP3  =  TENPS'*3^RTRN  '  1 .  •  ■  S  IGT^SC 
GO  TO  75 

tenfs,-tenpS'*!-::^rtrn  .  I .  >  ♦-): 

SB  '  IS;'  =TEMP3  +  SB  '  IS> 

COr^TINUE 

PPIr^T♦.■  NIT.  ELEMENT'S;'  TO  ELEMEflT  NIJMBEP='  »IS 

FPINT#,  ■CSB'KS>.KS  =  IS-'>  T-n  .  IS' 

iSDr<TINUE 

PPINT*, 

FPIr^T♦.■ 

PPINT*.' 

PETUPN 

Er^D 

SUBF  OUT  I  NE  STIFF  '  CR  •  S  PP  -  S  FZ  -  NZS  -  NFS  -  S  I G  T .  R  S  -  PNP  -  PNZ  -  LMR! '  •  NT"' 

D  I  ["'tEI I  ON  L  R  '  J  '  ''  *  J  M  ’  •  kP'  '  7  •  7  ^  E'  '  •  _  E’Z  '  7 ' 7  •  5  '  '  R  _  ■  5  ^  5  0.'  *  F  t'lF’  '  “c  i  7  <  ^  P 
riZ  '  ' 


THi:  FZL'INE  R'SEMBLE:  THE  FPCELEM  C  CEFF I S  I  Eri  r  DP  STIFFNEI!  MRTF  I 

DO  £0  J=1.NT 
DO  £i’i  f  =1  .  NT 
R : j '  f  '  =  ij .  i;i 
f  j=':' 

Tz  = 

-  'I  ■:  r  I  .  ft  - 

■  •  >  H  ’  ’ 


11 


ELSE  IF  '  H.  EC- 1 1  ■'  THEM 
NP  =  £' 

riz=i 

ME=1 

bO  TO  SO 

ELSE  IF  ai-EC-lE)  THEN 

NF=4 

HZ=£ 

GD  TO  SO 

ELSE  IF  '  H.  EO,  is:-'  THEtl 
NF-£ 

HZ=1 

ME=1 

GD  TO  £5 

ELSE  IF  ';:M.  EO.  14;:'  THEN 

np=s 

OZ^I 

60  TO  SO 

ELSE  IF  ':M.E0.15;>  THEM 

OF-S 

MZ=£ 

0E=1 

GO  TO  SO 

EL  IE  IF  'h'.EO.lG'  THEM 

r<p=4 

ffZ=£ 

GO  TO  cl 

ELSE  IF  'SLEO.l?'  THEM 
OF  =5 
riZ  =  £ 

NE=1 

GD  TO  S'S 

ELSE  IF  i:fLE0.1S;'  THEH 
riP=G 

riz=  ? 

GD  TO  ?0 


T  Cf  -  1  c  I  * *  M  r  f-  r '  '  !  .  '  *  ■  I  S  4 
i  M  •  ^  L'  '  *  ■?  *p  TOr^  '  1  .  -  I  tf 'F'C 


•  .  -  /-  T  .  .  D*-  .  I  Z  .  r  T 

■  V  -  .w7  '  . :  'P  ‘  .  GT .  '  ■•'u 

-  IT  ,  '  .  G'"' .  1  .  ■rpEO 

L  1  =4 


CM  L  L  r 


-C  ,  ,NC  .  *  :  F  z  '  I  Z  ■ 


b._ 


-CM- ■ cc .  jF-  OF' ♦IFZ ’  IZ'  J-*  O- ' 

GO  TO  SZ  _ 

TEM=:ff  'IF'  JF«  OP'  ♦-FZ  ■  ' 

b  0  T  0  c  , 

TCM=-PF'  IF,  IP.OF'^IFZ'  -F  ' 

TEOF  1  :.TEO^CG  •  JFI,  IH'  0  ■  ♦  ■  - 1  ■  -^OE^TEMP  1 


-c-.<Pt  ^TFbPi  »C-*H  'pi'  ■  1 

-  p'  -.iC  S  Tt  '  P  C'  <  _■  •  I  I'  ^  ^  ^  ' 

-  C  '  ’C'  ,  J  C  n  r-  t  ^  \ 


; ,  a  .  ♦»  -I 
■ ,  4  .  ♦!_  H 


111 


-HP:" 


el:e  if  'i_i'.Eij.-r' 
LE=:' 

1  Hr,  n 

C  ^  ^ 

EL  IE  'Lr.EO.-E' 

.  THEM 

L£  =  E 

ED  TC  -E 

EDDI- 

el:e  'JD.eC'.E' 

L 1  = 

IF  'Lr'.E'I'.Ci'  THEr4 

THE>^ 

ED  TC  IE 

EL  IE  IE  'LL.EO.l' 

THEM 

LE  =  4 

ED  TD  I- 

EL  IE  IF  'LD.EE.-I' 
LE  =  E 

■  THEM 

EC  TD  :--T 

END  IE 

EL  IE  !=■  'JZ.EE.::  ' 
L1=E 

IF  'LL. EE. O'  then 

THEtJ 

LE'=E 

EC  TC  :e 

el:e  if  'Lii.E!:'.;’ 

Ti-iFM 

LE=  ■ 

I'O  TC  '■? 

EL-E  IE  'LD.EE.i' 
LE  =  4 

edlif 

T  he  '  ■ 

Ef'DIF 

FE  Z  I  =  I  F'F  1  1  .  L  1  >  FfL. 

' ,  _i2_  . .  Pfc  ■  j; 

I  LE  '  FT'Z  '  IZ-Z  '  '  FHZ  ' 

I  Z-Z  :■  FMZ  '  I Z- 

ED  TC  4n 

1  i'  =  rC  I  -  JZ 

■ ’I  ~LZ "  ~  I Z 

rr  .  1  I'.  I'f  .  CD. 

,-T_  t  ,  :-n  Tq  a 

-  ■  ■' 

,  ,  ^  r  ~ 

I'-fZ  ■  ^7  r '  ■  r 

Z"  TC  F; 

f  t'  =  f<F  I-OF 

r.-j;,-r,c :  - :  p 
'F  ■  ^  Z.  PT.  [JC  ,  MI' 

.GT.d-Gu  TG  ^ 

njBPnijTir-iE  EGi'r^D  •  j':*  *  temp«h;::  ■ 

DiMErr:  rcr^  fp:  :  ■  -p:  r  - 


♦  TWT:  FOUTUiE  '.ELECT-;  Hhli  compute;  the  VhLUE;  OF  THE  EICUEIC 

♦  :  PL  I  ME  FUr'CTIOtJ;  hT  the  DUTEP  PHD  I  lit  hMD  top  'UPFhCE  of  the 

♦  F'^'DILEM  CL'LIMTEP.  THE-E  VRLIjET  HPE  PEOUIPEr  IM  THE  E  VELUPlT  I  DM  E- 

♦  THE  -I.iP'FHLE  IMTEhP'HL  TEP'M>-. 


L  i^fc : :  -  j: : 

LT^  j:  :-i:  ■ 

IF  'LI.EC.h.  THEfi 
Ll  =  l 

IF  'Lri.EC.O'  THEri 
LE=1 

CD  TO  40 

El:E  IF  -LD.EO.l-  TI-E.M 
CD  TO  4 1 j 

el:e  if  'LI'.ec.E'  theii 

LE=; 

CD  TO  40 
EMTIF 

EL  IE  IF  'Ll. EC.  I'  them 
L1=E 

IF  'LT.  THEM 


CO  TO  -0 

EL-E  IF 'Ll. EC.-!  ■ 

THEM 

LE=1 

CO  TO  4  1' 

EL  IE  I!^  '  LT.  EC  .  1  ■ 

THEM 

LE=; 

GO  TO  4  11 

Ef'DIF 

EL:e  IF  'LI.EC.E' 

THE!' 

Ll  =  ;' 

IF  'Lr'.EC.O'  THEM 

LE=,£ 

GO  TO  40 

el:e  if  ■Lr'.EC.-E'  theij 
LE-! 

EMTI'^ 

Ct.'l'l  T  C 

,  [1  i'  ■  ,  ,  z  r.'  ,  ‘.t-f- 

■r  ,  -  ,  _  -  ■ ,  ■  •  ■  -  -  t .  ■  : "  -  ■  ,  '  -  ;  ,  f-'  ■  1  .  -  f  ■ '  ■  '  ^  ' 

P  F.  ^'JF  M 


PPIHT^.'  MhTPIX  check  FDP  LIHCDMhL  DDMIHhMCE. 
ppIr^T♦.■ 

ppIr^T♦,  poi.i  ruftcahRL  element  cum  uf 

PPINT^. 

DD  40  1=1 -NT 
TEMF-0. 0 
nn  CO  j^i-NT 
temp=temp-»-hbc  'hc  a .  j:-  ■ 

TEMP  =  TEMP-FiEC  '  PC  I  -  I  :• 

F'  F’  I N  T  ♦  -  -  I  -  ■  »  H  j  '  I  -  I  '  -  -  T  E  M  F' 

ppIr^T♦- 

CONTINUE 

PPINT^. 

PPira^-'  MflTPIK  chect  FGP  CYMMETPY.  ■ 

FPINT^. 

DO  pn  1  =  1 -fiT 
DO  PO  J=l-I 

F'F'INT^-  Hj*-'  -  I-  -  -J-  '= 

1  C  •:  J  -  I  ' 

COr^TIr^LlE 

PETUPr^ 

END 

C  Ij  t;  F'  D IJ  T  I N  E  :  D  L  '■/  E  ■  1,1 ,  E ,  H  -  -  N  T  - :  P  ■' 

D I  r’Ert :  I Gr<  r  :  ■  c  o -  ■?  o  ■  - :  e  ■  ly '  - :  h  •  g  o  ■  - :  e  ■  5  o  • 


CUM  OF  PON  element:. 


-  _i  -  -  - 1  - 


•  THi:  PouTirtE  iclve:  the  ppdelem  mrtpice:  ec  the  itefrtt'e  met 

*  OF  lUCCEIdVE  DVEP-F'ELR:  RTIGN. 

FPIMT*, 

PPINT*, 

DO  CO  1  =  1- rd 
:  :r  '  I  '  =n. 

YE:  ■  I  ■  =  r'. 

IT=0 
IT=IT-k1 
DO  C.'i  ^=l-NT 
TEMF- :P ' K' 

-c  '  =  ‘  •  fd 

--■,‘0=T'rMp_c.-  ,  f  .  ;  . 

^  I  -  rp'^F  *i,'  ■  H  :  '  ^  -  f  '  '  L  '  r  ' 

DO  fr-o  1  =  1  -NT 
I P  '  'E'  '  I  '  .  El? .  0 .  CD  TO  P  0 
TDL=  ' '  p  ■  I  '  -‘<H '  I '  ■  :  :E:  •  1  • 

IF  '  HP  :  i  TGL  '  .  CT .  .  riO  1  ’  CO  TO  TO 
roriTiriijE 

Itp  jr’  ■:,ri 


rn  -:ii  [  =  i.r<T 


"CT*  ».  ’OLl'T’’Cr-  :  0- I-:- '--'E  D  .  THE  ■‘-I’'''--'.. 

tFir'T*-  THE  rKJMPEP  CP  ITEPCTIOrr:  I.IEPE  -it 
d^U'IFi'  ■■FC"’'^.C!T,c'j“r';;f,  iIOOFFIlIC 


PRINTS.  THE  PPOELEM  HhT  HCT  CONVERGEri  IN  '-IT.  ITEFhTIDN:. 
PRINT*. 

PPIr^T♦.  THE  EPPOP  =  .TDL 
PPINT*. 

PR  INT*.  '  .  NT  •  I  ■  .  1  =  1  .NT;:' 

RETURN 

END 

rUEPDUTINE  TFLUU  PNP .  PNZ .  XE .  P -  Z •  KT .  NZ3  .  NPS  .  TPHI  .  DPHI  .  Hp.  I  IGT.  HIH- 
ID  ■ 

DINEN'IDN  PNZ'-Z:Z:'  .PNP'-£:7-'  .’'P-GO:'  .  IZPF'4.4:' 

*  THi:  PDUTINE  CCMPUTEZ  THE  'ChTTEPED  hND  TOTRL  PERCTICN  PRTE 

*  FLUENCE. 

C  RLL  I ENRT  I Z  PF .  PNP .  FNZ .  R  .  Z .  NPZ  .  r^Z  Z  :■ 
f  Z  =  5 
TEMP  1  =  0. 

DC  4  i:i  1  =  1.4 
MZ=I IPF ■  I .  1  ■ 

IF  '  NZ,  GT.  rc:  'GG  TD  4  0 
►  Z  =  I  Z-1 

^  Rnf. 

Du  Z:  0  _i  =  1 , 4 
MC'=  I  :pp .  .1, 4 . 

IF'NP.  GT.rtPI  'GC  TO  jO 

riF  =  1  T*  .  '  11'='- 1  '  +r.p : *  •  nz- i  ■  ■ 

f  R=^  R-l 

TCMP  =  :  ,  rJF  '  *':PF  .  1  .  f  z.  PNZ  '  NZ-?''  .  PNZ  '  MZ-Z  '  .  PNZ  '  NZ-  1  '  .  PNZ  '  NZ  '  >  Z  ■  *  Z  ' 

1  1  .  ^  R  .  R' f R  *  f'' R  —  '  .  R  T'^R'  '  f‘1R'  — c  '  .  R’NR  '  f’iR'—  ^  '  .  P  f^R’  '  NR’  '  .  R  ' 

TEf’-Pl=TENPI  +TENF 
'ZCNT  irOJE 
iZCr^TirPJE 

H=  Z  C'PT  '  I  >:■  .  *HTHr‘  '  1  .  >  ' 

TPhI=m*TFmf1 

2H=Z-HP 

PZH=ZORT  'P**?:*ZH**Z  ' 

I  F  '  R  ZH  .  E L' .  * '  '  UC  TD  *r'  0 


-cl  ;  =  T  :  -  T  »  I  P  '  p  i:,  '  lj  .  •  -p-  ,  •  p,  ,  ,  -  p;  i  :  r, 

GC  TG  -O 

TRij=:  igt*e:  p  ' -z  hzH'*p 

IiFHI  ='r  D*E:  :P  '  -Thu  '  ■  lP*HTHr< .  i  .  .  ♦p-H**Z:' 

GG  TG 

pciriT*.  EPRCP.THE  Pul.  I'RUr^GT  PE  COr’IPUTED  RT  THE  P'.IPZT  PCJUT. 

-  Tr-cp. 

c  - '  p  f  I 

.  r, 

'  ,  p  c  PI  IT  r  *jir  T  cr-p. T  '  '  P  ‘  "C  .  PIJ  -  .  F  .  7  .  r-P  ^  ■  ■ 


HI  i  -p  ■;  UP  --E  -  I'U.’f  I'Z 


*  .'ELE'.  TED 


pciN-^ : . 


Du  4i.i  I  =  M,n^ 

IF'  Z.GE.PNZ' I>  .HhD.Z.LT.PMZ'  I  +  l  -'  .-bD  TU  5  0 
CUNT  I HUE 
DU  60  J=li4 
I  ZPF  ''J,  i  ;:.  =I4.J 
DU  70  I=0'HPZ-3 

I F  '  P .  GE  .  PHP  •:  I  >  .  Rr^D .  P .  LT .  PHP  a  +  1  >  :■  GU  TU  S  0 

CUHTIHUE 

DU  S5  J=1.4 

r:pp  ■  J,  =I+.J 

GU  TU  100 

PPIHT^.'  THE  ZPRlIRL  PUZITIUH  UF  -Z-  '  UP  -P- 
iPPUELEM  dimeh:  lurr: . 

'TUP 

PETUPN 


r:  HUT  I.IITHIH  THE 


Appendix  K 


Numerical  Results 

The  computer  subroutines  which  are  listed  in  Appendix.  J 
were  used  to  produce  numerical  results  for  varying  problem 
spatial  mesh  sizes  and  degrees  of  the  spherical  harmonic 
trial  function  expansion.  Some  of  these  results  are  presented 
in  Figures  10  through  44.  They  are  valid  for  the  problem  para¬ 
meters  which  were  presented  in  Chapter  V,  However,  the  cross- 
sections  which  were  used  do  not  accurately  represent  the 
values  for  air  at  sea  level.  Also,  the  air-ground  interface 
was  not  included  in  the  problem  domain  and  therefore  all 
ground  effects  xvere  ignored. 

Decausc  of  the  time  constraints  on  this  research  project 
neither  an  evaluation  of  the  accuracy  of  these  results  nor 
a  comparison  to  a  discrete  ordinate  or  Monte  Carlo  calcula¬ 
tion  was  accomplished.  Therefore,  the  results  are  presented 
solely  in  an  attempt  to  show  that  finite  element  space- 
angle  synthesis  is  a  viable  solution  technique  for  solving 
the  two-dimensional  steady  state  anisotropic  Boltzmann 
equation.  They  are  not  meant  to  represent  a  precise  and 
exact  solution  to  the  a i r-ovc r-ground  problem,  but  rather 
to  demonstrate  that  FF.SAS  may  be  a  feasible  alternate  solution 
technique  to  Monte  Carlo  and  discrete  ordinates. 
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